Illumina pair end read extraction.

[JDinis]

Hi all,

I am working with a data set of 161x161 pair-end sequences generated on the Illumina MiSeq. In short, my goal is to extract all pair-end reads that overlap a 500 bp window of my assembly (easy part). I want to be able to show linkages of polymorphisms within this window.

The problem that I am having is that I also want to extract all paired reads that are associated to different polymorphisms, such as all reads that overlap positions 700, 750 and 1005 based on my MSA (hard part).

Is there a good program that can extract reads from MSA, SAM or BAM alignments based on the qualifier of nucleotide positions? If so what are your suggestions?

Thanks all for your time.

JD

[SeekAnswers]

You can specify the coordinates that you like to extract from your bam with samtools view and make a smaller bam file.

Then use something like bam2fastq to pull out the reads in that region.

[JDinis]

To my knowledge of samtools (correct me if I am wrong) that function only work when extracting a range, extract all reads from 400 to 900 nt based on a reference. This will not work for my application.

My goal is to extract paired reads that overlap coordinates 700 and 750 and 1000, it is similar in concept but a very different output. However, if I am wrong could you please let me know and could you please explain the syntax for such a search in samtools.

Thanks for the reply.

[Dario1984]

Use the scanBam function from the R package Rsamtools. You can give it a set of ranges of width 1, using the what parameter to ScanBamParam.