sWGS检测CNV的一点探索

2021-04-19 14:57:33 浏览数 (2)

试用了两个软件用于测试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这种测序深度。

na panel

0 人点赞