This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Paired-end reads that span a position in an assembly

I need to find pairs of reads that span a certain position in an assembly. However, neither needs to necessarily correspond to that position. To be clear, here is an example:

assembly: ----------------------[=====================]----------------------------

pair A:   -----[<<<<<]------------------------------------[>>>>>]------------------

pair B:   -------------------------[<<<<<]------------------------[>>>>>]----------

I am interested in the positions [=====] in the assembly. I am looking for pairs of reads like read pair A. I can find reads like read pair B using bedtools intersect, but I am not interested in these.

bedtools genome assembly paired-end
• 1,159 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by saladi • saladi

1 answer

using samjdk. http://lindenb.github.io/jvarkit/SamJdk.html Say you want spanning RF01:1000-2000.

samtools view -u S1.bam RF01 | java -jar dist/samjdk.jar -e 'if(record.getReadUnmappedFlag() || !record.getReadPairedFlag() || record.getMateUnmappedFlag() || !record.getContig().equals(record.getMateReferenceName()) || !record.getContig().equals("RF01")) return false;  return(Math.min(record.getStart(),record.getMateAlignmentStart()) < 1000 && Math.max(record.getStart(),record.getMateAlignmentStart())>2000) ;'
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this