Hot answers tagged

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: ...
Devon Ryan's user avatar
  • 19.6k
6 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 ...
bli's user avatar
  • 3,130
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 ...
Maximilian Press's user avatar
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 =...
Kyle's user avatar
  • 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 ...
Devon Ryan's user avatar
  • 19.6k
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 ...
John Marshall's user avatar
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) ...
terdon's user avatar
  • 9,846
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: ...
Devon Ryan's user avatar
  • 19.6k
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 ...
Devon Ryan's user avatar
  • 19.6k
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 ...
Devon Ryan's user avatar
  • 19.6k
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: ...
Devon Ryan's user avatar
  • 19.6k
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
The Unfun Cat's user avatar
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') ...
Bioathlete's user avatar
  • 2,574
3 votes
Accepted

Is it possible to, all in memory without writing a file, given a list of AlignedSegments and a bam header, build a bam, index and get coverage?

Assuming you are on an operating system that supports them, like Linux, you can try using named pipes in a ramdisk or tmpfs. Not the most elegant approach, but I don't know if what you are asking for ...
terdon's user avatar
  • 9,846
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 ...
RHS's user avatar
  • 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 ...
Ian Sudbery's user avatar
  • 3,321
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 ...
finswimmer's user avatar
  • 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-...
Karel Břinda's user avatar
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 ...
Bioathlete's user avatar
  • 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 ...
opplatek's user avatar
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 ...
StupidWolf's user avatar
  • 1,688
2 votes
Accepted

How to subsample BAM based on read line number?

Samtools+awk: ...
user172818's user avatar
  • 6,525
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 ...
winni2k's user avatar
  • 2,256
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 ...
bricoletc's user avatar
  • 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 ...
Mehrdad Bakhtiari's user avatar
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 ...
Daniel Standage's user avatar
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,...
gringer's user avatar
  • 13.9k
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: ...
Scott Gigante's user avatar

Only top scored, non community-wiki answers of a minimum length are eligible