9
votes
Accepted
Quantifying reads mapping to multiple loci
You almost had the correct python code already, you just need to filter out secondary alignments:
...
- 19.5k
5
votes
Is there an efficient way to extract CIGAR strings for read pairs from bam files with python?
I think that the easiest way is to work with query-name sorted groups of reads. In that case, mates will be adjacent in the sort and you can use that to extract the paired CIGARs.
If you are depending ...
- 3,649
5
votes
Accepted
How to filter a SAM file by a bed file?
Thanks to the helpful comments I figured it out:
samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam
-b = Output as bam file
-h = include header
-L =...
- 171
5
votes
pysam or piping samtools view to a python script
It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will ...
- 19.5k
5
votes
Accepted
What does "fetching by region is not available for SAM files" mean?
Your 3_Tms_1_mapped.bam file, despite its filename extension, is in fact a bgzipped SAM file. You can verify this using htsfile, which is a small utility packaged ...
- 946
5
votes
What does "fetching by region is not available for SAM files" mean?
That isn't actually a bam file as John Marshall figured out. I am keeping the rest of my answer since it could be useful to someone else, but the issue here was that you had a compressed (bgzipped) ...
- 9,452
5
votes
Accepted
Using pysam with cython: htslib/kstring.h not found
Based on a suggestion by Devon Ryan, I searched how to set C_INCLUDE_PATH for cython, and found this issue. It refers to a (currently not working) feature present ...
- 3,090
4
votes
Accepted
Filtering bases based on phred qualities with pysam
The PileupRead object as a query_position attribute, which you can use for this:
...
- 19.5k
4
votes
Accepted
Extracting the CIGAR string from a BAM via Python?
If you really do just want the cigar string then it's read.cigarstring. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless ...
- 19.5k
4
votes
Access base aligned to particular reference position
You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased ...
- 19.5k
4
votes
Accepted
How to convert BAM file to bigWig in Python?
deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:
...
- 19.5k
3
votes
How to convert BAM file to bigWig in Python?
pyranges:
import pyranges as pr
gr = pr.read_bam("your_file.bam")
gr.to_bigwig("out.bw", chromosome_sizes=pr.data.chromsizes()) # for hg19
- 451
3
votes
filter secondary alignments using pysam
Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the read.has_tag('XS') ...
- 2,574
3
votes
Access base aligned to particular reference position
I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)
This code takes advantage of the cigar string. It returns None if the base has been ...
- 131
2
votes
Chunk alignment in a name sorted bam for parallel processing
I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing ...
- 3,271
2
votes
Chunk alignment in a name sorted bam for parallel processing
python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is ...
- 1,342
2
votes
How do I rewrite a read group using pysam?
If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-...
- 1,899
2
votes
Accepted
pysam pileup: what reads appear in the pileup?
The two entries in the sam file represent mate pairs for the given ID. You can tell the difference based on the sam flags. Picard has a nice tool to determine the meaning of the flags.
You are ...
- 2,574
2
votes
Convert paired-end BAM into a single-end BAM and keep all the reads
MRNM stands for "Mate reference index". So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.
The ...
- 21
2
votes
Convert paired-end BAM into a single-end BAM and keep all the reads
Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.
You can also use samtools rmdup too
...
- 1,678
2
votes
Accepted
2
votes
Accepted
Is there a tool that can perform a read-group-aware mpileup from a single file?
I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming_pileupy and it's available ...
- 2,216
2
votes
How to filter a SAM file by a bed file?
You can also do this with bedtools intersect:
bedtools intersect -abam input.bam -b bedfile.bed > output.bam
This should be ...
- 161
1
vote
Accepted
Find all the bases for given reference position
pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups ...
1
vote
Aligned base strand in pysam `pileupcolumn`
As @terdon mentions, it's certainly possible to have an A aligned on both strands. If a read has a mismatch at the location of interest and that mismatch is the ...
- 5,060
1
vote
Access base aligned to particular reference position
This is not [yet] a complete answer, but hopefully it helps get you on the right path.
I wrote some code to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference,...
- 13.5k
1
vote
Reject reads with low quality bases from a Bam file through pysam
For can do this by accessing the basecall qualities from PileupRead.alignment. For example:
...
- 2,133
Only top scored, non community-wiki answers of a minimum length are eligible