前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >sWGS检测CNV的一点探索

sWGS检测CNV的一点探索

作者头像
用户1075469
发布2021-04-19 14:57:33
1.3K0
发布2021-04-19 14:57:33
举报
文章被收录于专栏:科技记者

试用了两个软件用于测试CNV检测,虽然没有取得什么结果,记录和分享一下。

ichorCNA笔记

这个软件可以检测切除的肿瘤组织,识别其中的肿瘤细胞含量,也可以用来检测纯肿瘤组织。可以有参考,也可以不用,官方提供了参考,可以自建。

1、 软件安装

软件官网:https://github.com/broadinstitute/ichorCNA

代码语言:javascript
复制
library(devtools)
install_github("broadinstitute/ichorCNA")

2、软件使用

代码语言:javascript
复制
# 1、准备数据,分块10K
hmmcopy_utils/bin/readCounter --window 10000 --quality 20 \
    --chromosome "chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX,chrY" \
/media/vd1/fastq/WGS/WGS-CNV/fastq/75bp_trim/sLK-1054_75_sorted.bam > 1054_75.wig
# 2、应用
Rscript ./runichroCNA.R --id 1054_75 \
 --WIG 1054_75.wig --ploidy "c(2)" --normal "c(0.5,0.6,0.7,0.8,0.9)" --maxCN 5 \
  --gcWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/gc_hg19_10kb.wig \
  --mapWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/map_hg19_10kb.wig \
  --centromere /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt \
  --normalPanel /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/HD_ULP_PoN_1Mb_median_normAutosome_mapScoreFiltered_median.rds \
  --includeHOMD False --chrs "c(1:22)" --chrTrain "c(1:22)" \
  --estimateNormal True --estimatePloidy True --estimateScPrevalence True --txnE 0.9999 --txnStrength 10000 \
  --scStates "c()"  --outDir ./

参数说明:

1) --id 样本名

2)--WIG 前一步准备的样本文件

3) --gcWig GC含量的文件,illumina测序对GC含量有偏好,需要去除影响

4) --mapWig 这是参考基因组的10K window图谱

5) --centromere 每个染色体都有的着丝粒,去除

6) --normalPanel 阴性参考,可选,有的话去噪效果好,可自己用阴性构建

7) --includeHOMD 特别大Bin时如1M需要

8) --estimateScPrevalence FALSE --scStates "c()" 针对体细胸样本时需要,亚克隆 9) --chrs "c(1:22)" --chrTrain "c(1:22)" 只训练和分析常染色体

10) --estimateNormal True 评估正常的,这个参数默认,可以不加

11) --estimatePloidy True 评估肿瘤细胞倍性,也是默认,可不加

12) --estimateScPrevalence True 亚克隆情况,体细胞需要,可设置False

3、结果

以3个阴性样本为对照,生成Normal 参考,去噪,衡量两个阳性样本,发现只有一个缺失较大片段(20K+,s1054)的结果可以看出,只是图不太一样,并不能看到明确缺失。

4、参考Panel构建

可以使用阴性对照建立自己的基线参考,不限样本,越多越好,去噪有一定效果,这里采用10K窗口进行构建:

代码语言:javascript
复制
#panel
Rscript createPanelOfNormals.R \
    --filelist ./noraml.txt \
    --gcWig  /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/gc_hg19_10kb.wig \
    --mapWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/map_hg19_10kb.wig \
    --centromere /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt \
    --outfile normal_sy

二、 QNDASeq笔记

1、软件安装

代码语言:javascript
复制
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("QDNAseq")
BiocManager::install("QDNAseq.hg19")

2、软件使用

代码语言:javascript
复制
#增加内存限制,防止不足,不设置会报错
options(future.globals.maxSize=5000000000)
library(QDNAseq)
bins <- getBinAnnotations(binSize=5)#获取每一个bin的基因注释,联网下载
bins
readCounts <- binReadCounts(bins,bamfile='sLK-1054_75_sorted.bam')
readCounts
pdf("readCounts.pdf")
plot(readCounts, logTransform=FALSE, ylim=c(-50, 200))
highlightFilters(readCounts, logTransform=FALSE,residual=TRUE, blacklist=TRUE)
readCountsFiltered <- applyFilters(readCounts,residual=TRUE, blacklist=TRUE,chromosomes=c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","18","19","20","21","22","X", "Y"))
dev.off()
pdf("mapping.pdf")
isobarPlot(readCountsFiltered)
readCountsFiltered <- estimateCorrection(readCountsFiltered, ncpus=8)
dev.off()
pdf("sd.pdf")
noisePlot(readCountsFiltered)
copyNumbers <- correctBins(readCountsFiltered)
copyNumbers
copyNumbersNormalized <- normalizeBins(copyNumbers)
copyNumbersSmooth <- smoothOutlierBins(copyNumbersNormalized)
dev.off()
pdf("cnr.pdf")
plot(copyNumbersSmooth)
exportBins(copyNumbersSmooth, file="sample.txt")
exportBins(copyNumbersSmooth, file="sample.igv", format="igv")
exportBins(copyNumbersSmooth, file="sample.bed", format="bed")
dev.off()
######1.3Downstream analyses#########
copyNumbersSegmented <- segmentBins(copyNumbersSmooth, transformFun="sqrt")
copyNumbersSegmented <- normalizeSegmentedBins(copyNumbersSegmented)
pdf("cns.pdf")
plot(copyNumbersSegmented)
copyNumbersCalled <- callBins(copyNumbersSegmented)
dev.off()
pdf("cnc.pdf")
plot(copyNumbersCalled@featureData@data)
exportBins(copyNumbersCalled, format="vcf")
exportBins(copyNumbersCalled, format="seg")
dev.off()

总结,没有可靠结果,软件暂时认为不可用于10Kb以下的CNV检测,特别是1x这种测序深度。

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-03-28,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 科技记者 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 试用了两个软件用于测试CNV检测,虽然没有取得什么结果,记录和分享一下。
  • ichorCNA笔记
    • 1、 软件安装
      • 2、软件使用
        • 3、结果
          • 4、参考Panel构建
          • 二、 QNDASeq笔记
            • 1、软件安装
              • 2、软件使用
              领券
              问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档