Question: calculate fragments depth with pysam
1
gravatar for electronsource75
2.5 years ago by
electronsource7510 wrote:

Hello, I am currently working with pysam to calculate sequencing “depth” of fragments. I use a bam file with Pair End reads and I managed to create a function that extracts the fragment name, length and the start and stop. I use python and the pysam for this project

My problem lies in two points :

the first is to find a better way to calculate and store the fragments. Because the most efficient way is to use a bam file ordered by read names instead of genomic coordinates that way the read is always followed by its mate. But in pysam, I cant use a bam file sorted by read names because it is impossible to create an index (with samtools) for this kind of file and pysam need an index

the second point is to use the data stored to compute the fragment's depth. I ’ll calculate the fragment's depth for each chromosomes using multi-threading.

pysam samtools python genome • 1.7k views
ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by electronsource7510

why do you need the mate ?

ADD REPLYlink written 2.5 years ago by Pierre Lindenbaum131k

because with the mate I can have the position where the fragment ends

ADD REPLYlink modified 2.5 years ago • written 2.5 years ago by electronsource7510

isn't RPOS (mate pos) enough ? sometimes one can also find the mate cigar string in the read metadata.

ADD REPLYlink written 2.5 years ago by Pierre Lindenbaum131k

Thank you I search through pysam and found an option to get the rpos

ADD REPLYlink written 2.5 years ago by electronsource7510

bamCoverage from deepTools would directly do this for you.

ADD REPLYlink written 2.5 years ago by Devon Ryan97k

yes it does the job pretty well but I want if possible do my own program using pysam to do the same thing

ADD REPLYlink written 2.5 years ago by electronsource7510

Hi, I a have another question when I observed the fragments length some of them have a negative value like (-54). Why? I use pysam and I only took reads which are paired and correctly mapped.

ADD REPLYlink written 2.5 years ago by electronsource7510
1

Please see the SAM specification.

ADD REPLYlink written 2.5 years ago by Devon Ryan97k
1
gravatar for Devon Ryan
2.5 years ago by
Devon Ryan97k
Freiburg, Germany
Devon Ryan97k wrote:

For recent versions of pysam and assuming you have PE reads it's simple enough to just use read.template_length to get the full fragment length. If you need to handle SE reads or older pysam versions then you can get slightly more elaborate as shown here. You'll then want to ignore one of the mates (unless they're not properly paired, which you may want to treat specially (in deepTools we would use the median population fragment length in that case)).

ADD COMMENTlink written 2.5 years ago by Devon Ryan97k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1101 users visited in the last hour