sWGS检测CNV的一点探索

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

ichorCNA笔记

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

1、 软件安装

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

  1. library(devtools)

  2. install_github("broadinstitute/ichorCNA")

2、软件使用

  1. # 1、准备数据,分块10K

  2. hmmcopy_utils/bin/readCounter --window 10000 --quality 20

  3. --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"

  4. /media/vd1/fastq/WGS/WGS-CNV/fastq/75bp_trim/sLK-1054_75_sorted.bam > 1054_75.wig

  5. # 2、应用

  6. Rscript ./runichroCNA.R --id 1054_75

  7. --WIG 1054_75.wig --ploidy "c(2)" --normal "c(0.5,0.6,0.7,0.8,0.9)" --maxCN 5

  8. --gcWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/gc_hg19_10kb.wig

  9. --mapWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/map_hg19_10kb.wig

  10. --centromere /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt

  11. --normalPanel /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/HD_ULP_PoN_1Mb_median_normAutosome_mapScoreFiltered_median.rds

  12. --includeHOMD False --chrs "c(1:22)" --chrTrain "c(1:22)"

  13. --estimateNormal True --estimatePloidy True --estimateScPrevalence True --txnE 0.9999 --txnStrength 10000

  14. --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窗口进行构建:

  1. #panel

  2. Rscript createPanelOfNormals.R

  3. --filelist ./noraml.txt

  4. --gcWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/gc_hg19_10kb.wig

  5. --mapWig /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/map_hg19_10kb.wig

  6. --centromere /media/vd1/zjd/Miniconda3/lib/R/library/ichorCNA/extdata/GRCh37.p13_centromere_UCSC-gapTable.txt

  7. --outfile normal_sy

二、 QNDASeq笔记

1、软件安装

  1. if (!requireNamespace("BiocManager", quietly = TRUE))

  2. install.packages("BiocManager")

  3. BiocManager::install("QDNAseq")

  4. BiocManager::install("QDNAseq.hg19")

2、软件使用

  1. #增加内存限制,防止不足,不设置会报错

  2. options(future.globals.maxSize=5000000000)

  3. library(QDNAseq)

  4. bins <- getBinAnnotations(binSize=5)#获取每一个bin的基因注释,联网下载

  5. bins

  6. readCounts <- binReadCounts(bins,bamfile='sLK-1054_75_sorted.bam')

  7. readCounts

  8. pdf("readCounts.pdf")

  9. plot(readCounts, logTransform=FALSE, ylim=c(-50, 200))

  10. highlightFilters(readCounts, logTransform=FALSE,residual=TRUE, blacklist=TRUE)

  11. 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"))

  12. dev.off()

  13. pdf("mapping.pdf")

  14. isobarPlot(readCountsFiltered)

  15. readCountsFiltered <- estimateCorrection(readCountsFiltered, ncpus=8)

  16. dev.off()

  17. pdf("sd.pdf")

  18. noisePlot(readCountsFiltered)

  19. copyNumbers <- correctBins(readCountsFiltered)

  20. copyNumbers

  21. copyNumbersNormalized <- normalizeBins(copyNumbers)

  22. copyNumbersSmooth <- smoothOutlierBins(copyNumbersNormalized)

  23. dev.off()

  24. pdf("cnr.pdf")

  25. plot(copyNumbersSmooth)

  26. exportBins(copyNumbersSmooth, file="sample.txt")

  27. exportBins(copyNumbersSmooth, file="sample.igv", format="igv")

  28. exportBins(copyNumbersSmooth, file="sample.bed", format="bed")

  29. dev.off()

  30. ######1.3Downstream analyses#########

  31. copyNumbersSegmented <- segmentBins(copyNumbersSmooth, transformFun="sqrt")

  32. copyNumbersSegmented <- normalizeSegmentedBins(copyNumbersSegmented)

  33. pdf("cns.pdf")

  34. plot(copyNumbersSegmented)

  35. copyNumbersCalled <- callBins(copyNumbersSegmented)

  36. dev.off()

  37. pdf("cnc.pdf")

  38. plot(copyNumbersCalled@featureData@data)

  39. exportBins(copyNumbersCalled, format="vcf")

  40. exportBins(copyNumbersCalled, format="seg")

  41. dev.off()

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


本篇文章来源于微信公众号: 微因

发表评论