Hot answers tagged

9

You almost had the correct python code already, you just need to filter out secondary alignments: def get_reads_hist(bam): bam = pysam.AlignmentFile(bam, 'rb') counts = Counter() for query in bam.fetch(): if query.is_secondary: continue nh_count = Counter(dict(query.get_tags())['NH']) counts += nh_count ...


5

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 on position sorting to extract the alignments of interest efficiently, you can maybe do the following: extract alignments in the region of interest plus a lot ...


5

deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is: from deeptools import writeBedGraph from deeptools.getScaleFactor import get_scale_factor wr = writeBedGraph.WriteBedGraph(options and input files) func_args = {'scaleFactor': get_scale_factor(...options...)} wr.run(writeBedGraph.scaleCoverage, ...


5

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 with HTSlib: $ htsfile 3_Tms_1_mapped.bam 3_Tms_1_mapped.bam: SAM version 1.3 BGZF-compressed sequence data (For files that really are in BAM format, it reports BAM version 1 compressed sequence ...


5

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) sam file and not an actual bam file and that's why you were getting that error. When I sorted your file in preparation for indexing it, I converted to a bam which ...


5

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 not be the bottleneck in your processing, rather the python script itself will be.


5

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 in cython's documentation, that should let libraries inform cython of the correct localisation for the header files in the module's setup: ext_modules = cythonize( "bam25prime/libbam25prime....


4

The PileupRead object as a query_position attribute, which you can use for this: for pileupcolumn in samfile.pileup() : for pileupread in pileupcolumn.pileups: if (pileupread.alignment.mapping_quality <= 15): continue if not pileupread.is_del and not pileupread.is_refskip: if pileupread.alignment....


4

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 you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character ...


4

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 calls. for read in bam.fetch(some_position): positions = read.get_reference_positions(): if some_other_position in positions: idx = positions.index(some_position) idx2 = ...


3

pyranges: import pyranges as pr gr = pr.read_bam("your_file.bam") gr.to_bigwig("out.bw", chromosome_sizes=pr.data.chromsizes()) # for hg19


3

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') continue block. I think that if you filter out is_secondary and is_supplementary you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly. for read in bam....


3

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 deleted or clipped, or if the bam file did not store the sequence. from typing import Optional import pysam def find_base_in_alignment(alignment: pysam....


2

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 code if going the be the slow bit. import pysam infile = pysam.AlignmentFile("input.bam") chunk_size = 10000000 outfile_pattern = "output_segement%d.bam" chunk =...


2

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 reached the next chunk can only created of the read name is different. $ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" ' BEGIN { part=0; line=n } /^@/ {header = header$0"\n"; ...


2

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-overview.html#AddOrReplaceReadGroups. Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not ...


2

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 missing coverage within the first read. I think that it is related to the low quality bases. samtools mpileup filters bases with a quality below 13 and pysam just ...


2

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 problematic line in your code is: line.next_reference_id = 0 This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (next_reference_id). ...


2

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 samtools rmdup <bamfile> <output> Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads. I would suggest removing the duplicates and then ...


2

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 from pypi for installation with pip. After installation the command would be: spileup input.bam sample_names.txt The tool is pre-alpha and missing some features....


1

pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting ...


1

Thanks to the helpful comments I figured it out: samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam -b = Output as bed file -h = include header -L = only output alignments overlapping the input bed file


1

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 complement of the reference base, then you are likely to see A in both strands. However, it is very unlikely you will see the same base at high coverage in both strands. For example, if your total ...


1

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, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased): my $refPos = $F[3]; my $cigar = $...


1

For can do this by accessing the basecall qualities from PileupRead.alignment. For example: _BASE_QUAL_CUTOFF = 30 with pysam.Samfile(in_file, 'rb') as bam_in: for column in bam_in.pileup(): for read in column: qual = read.alignment.query_qualities if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]): ...


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