So currently I have been using cn.mops for CNV detection. The result for one sample is confusing since there is overlap between different CNV regions in cnvr output.
It reads as follows:
seqnames | start | end | width | strand | X1315_1.final.bam | |
1 | 1 | 144813741 | 249213345 | 104399605 | * | CN3 |
2 | 2 | 13140524 | 13150289 | 9766 | * | CN2 |
3 | 3 | 44143016 | 44170147 | 27132 | * | CN2 |
4 | 4 | 19526785 | 19545625 | 18841 | * | CN3 |
5 | 7 | 69864128 | 69924740 | 60613 | * | CN2 |
6 | 7 | 73050706 | 124472685 | 51421980 | * | CN4 |
7 | 7 | 86116372 | 115099423 | 28983052 | * | CN2 |
8 | 8 | 19904258 | 59107586 | 39203329 | * | CN2 |
It appeared to be an overlapping in chr7 with [86116372,115099423] assigned with CN2 while in the previous line it was given CN4. It's confusing, so I re-check the cnvs output.
It reads:
21 | 7 | 38217807 | 38937729 | 719923 | * | 1315_1.final.bam | -1.000003911 | -2.0535 | CN1 |
22 | 7 | 56496077 | 56949839 | 453763 | * | 1315_1.final.bam | -0.472927838 | -1.7989 | CN1 |
23 | 7 | 121946569 | 121950131 | 3563 | * | 1315_1.final.bam | 5.887672454 | 4.6284 | CN4 |
It's not overlapping however it doesn't exactly consistent with cnvr output.
SO the question is: which output should I use? Is there an explanation for the difference between cnvr and cnvs? what may be the cause of the cnvr results overlapping?
Here is the command I use:
library(cn.mops)
BAMFiles <- c('/home/export/data/P101SC17020317-01juyongzhi/Results/Bam/1315_1.final.bam','/home/export/data/P101SC17020317-01juyongzhi/Results/Bam/1315_2.final.bam','/home/export/data/P101SC17020317-01juyongzhi/Results/Bam/sq1315_pbmc.final.bam')
segments <- read.table('/home/export/data/1194/bam/result_bed_no_chr.txt',sep="\t",as.is=TRUE)
segments <- unique(segments)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr)
resRef <- referencecn.mops(cases=X[,1],controls=X[,3],classes=c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6","CN7","CN8","CN16","CN32","CN64","CN128"),I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64),segAlgorithm="DNAcopy")
resRef <- calcIntegerCopyNumbers(resRef)
segm <- as.data.frame(segmentation(resRef))
CNVs <- as.data.frame(cnvs(resRef))
CNVRegions <- as.data.frame(cnvr(resRef))
write.csv(segm,file="/home/export/CG/copynumber_clonal_evolution/result/1315/1315_1_segmentation.csv")
write.csv(CNVs,file="/home/export/CG/copynumber_clonal_evolution/result/1315/1315_1_cnvs.csv")
write.csv(CNVRegions,file="/home/export/CG/copynumber_clonal_evolution/result/1315/1315_1_cnvr.csv")
Thanks for the answer. However, I only use cn.mops for only one sample (1315_1). So there should not be a second individual.