Hello,
I am wondering why the MEDIPS.createSet() uses isFirstMateRead=T and isSecondMateRead=F when paired=TRUE?
For example, I have a PE bam file that contains only properly paired and mapped reads with 54,073,380 reads(~27M R1 and 27M R2).
When I import with:
sample <- MEDIPS.createSet(file="sample.bam",BSgenome="BSgenome.Hsapiens.UCSC.hg38",uniq=1e-3,shift=0,window_size=100,chr.select=c(paste0("chr",1:22,"chrX","chrY"),paired=T)
Reading bam alignment sample.bam
Selecting chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY
Total number of imported first mate reads in properly mapped pairs: 20772475
scanBamFlag: isPaired = T, isProperPair=TRUE , hasUnmappedMate=FALSE, isUnmappedQuery = F, isFirstMateRead = T, isSecondMateRead = F
Mean insertion size: 196.5556 nt
SD of the insertion size: 84.0891 nt
Max insertion size: 592 nt
Min insertion size: 51 nt
Creating GRange Object...
Keep at most 1 first mate read(s) mapping to the same genomic location
Number of remaining short reads: 18634486
Calculating genomic coordinates...
Creating Granges object for genome wide windows...
Calculating short read coverage at genome wide windows...
When I import with:
sample <- MEDIPS.createSet(file="sample.bam",BSgenome="BSgenome.Hsapiens.UCSC.hg38",uniq=1e-3,shift=0,window_size=100,chr.select=c(paste0("chr",1:22,"chrX","chrY"),paired=F)
Reading bam alignment sample.bam
Selecting chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY
Total number of imported short reads: 35457318
Creating GRange Object...
Keep at most 1 read(s) mapping to the same genomic location
Number of remaining reads: 29655862
Calculating genomic coordinates...
Creating Granges object for genome wide windows...
Calculating short read coverage at genome wide windows...
And I cannot change the settings within the function to import F and R reads:
sample <- MEDIPS.createSet(file="sample.bam",BSgenome="BSgenome.Hsapiens.UCSC.hg38",uniq=1e-3,shift=0,window_size=100,chr.select=c(paste0("chr",1:22,"chrX","chrY"),paired=T,isSecondMateRead=T)
Error in MEDIPS.createSet(file = hivfile, BSgenome = BSgenome, uniq = uniq, :
unused argument (isSecondMateRead = T)
Is there a biological reason why we would not want to use the R2 reads in MEDIPS? If not- how can I utilize both F and R reads in the analysis?
Thank you in advance, agm
