refSeqName issue in cn.MOPS
1
0
Entering edit mode
ashkot09 • 0
@ashkot09-9452
Last seen 5.7 years ago

I was wondering if someone could help.

I recently came across cn.MOPS and it's very good. While I was trying to use it I noticed that it asks for refSeqName and initially it led me to believe it is looking for the name of the reference file. But I don't think that is the case so then what should be used there:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/alignment/HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20120522.bam ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/alignment/HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20120522.bam.bai

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00099/alignment/HG00099.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00099/alignment/HG00099.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam.bai

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00100/alignment/HG00100.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00100/alignment/HG00100.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam.bai

and used to following code:

startMops <- function() {

library(cn.mops)

BAMFiles <- list.files(pattern=".bam\$")

res <- cn.mops(bamDataRanges)

plot(res,which=1)

data(cn.mops)

result <- calcIntegerCopyNumbers(cn.mops(XRanges))

segm <- as.data.frame(segmentation(result))

cnvs <- as.data.frame(cnvs(result))

cnvregions <- as.data.frame(cnvr(result))

write.csv(segm,file="segmentation.csv")

write.csv(cnvs,file="cnvs.csv")

write.csv(cnvregions,file="cnvr.csv")

return(1)

}

but i get error about the refSeqName.

Can you please let me know what I could be using.

CNV cn.mops • 1.1k views
2
Entering edit mode
@gunter-klambauer-5426
Last seen 8 months ago
Austria

Dear Ashwin,

The parameter "refSeqName"  relates to the names of the chromosomes as they are named in the BAM file header. Typically, these are "chr1", "chr2",... or "1", "2",... depending on the version of the reference genome.

If you want to analyze chromosome 1, then use:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,refSeqName="chr1",mode="unpaired",WL=5000)

If you want to analyze chromosome 1-22, X and Y, then use:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,refSeqName=c(paste0("chr",1:22),"chrX","chrY"),mode="unpaired",WL=5000)

Alternatively, you can try:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,refSeqName="1",mode="unpaired",WL=5000)

bamDataRanges <- getReadCountsFromBAM(BAMFiles,refSeqName=c(paste0("",1:22),"X","Y"),mode="unpaired",WL=5000)

Otherwise, your code seems to be okay. I would recommend to use more samples than just 3.

Regards,

Günter

0
Entering edit mode

And for what it's worth the sequence names ('levels') in the bam files can be discovered with seqlevels(BamFileList(BAMFiles)).

0
Entering edit mode

where is the method seqlevels located?

0
Entering edit mode

The function  "seqlevels" is a function of "GenomeInfoDb" and should be automatically available when you load cn.mops by "library(cn.mops)".

0
Entering edit mode

0
Entering edit mode

Thanks Gunter, this has worked.