15

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you'...


14

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.


13

Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and runs it robustly. Snakemake is my favourite. It's (nearly) pure Python and can ...


12

You can just set the fields you don't need to *: samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam - This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields. The resulting BAM ...


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 ...


9

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.


9

It is still slow but grep has a -f option to take in a file samtools view inbam.bam | grep -f read_names.txt > read_locs.txt


9

one-liner Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence. Note, this one-liner does not double-count reads because of uniq. samtools view file_sorted.bam | \ ...


8

Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency. So samtools merge (use the most up-to-date version, as there have been improvements in merge header handling in the 1.3.x and 1.4.x versions), picard MergeSamFiles, etc. These tools need to hold all the input BAM files ...


8

I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", ...


7

Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case). For use in featureCounts my advice is to get ...


7

I used this in the past for ChIP-seq data and it generated SNVs: samtools mpileup \ --uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \ --fasta-ref ref.fa file.bam \ | bcftools call --consensus-caller \ > out.vcf This was samtools 1.3 in case that makes a difference.


7

A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block so that the whole file does not need to be decompressed and therefore can be accessed randomly. This isn't true of files compressed as a single block using gzip....


6

I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, bioD, bioGo or rust-bio). Formatting an entire SAM is fairly expensive. You can see this by comparing samtools view aln.bam > /dev/null and samtools view -u aln....


6

The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power. As this question suggests, if ...


6

You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file. bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov" You can also provide a BED file if you just want to calculate in the target region.


6

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


6

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the --split RG argument. Here's some example output from a BAM file with combined read groups: $ samtools stats --...


6

Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by making a small BAM file either with or without duplicate entries but some marked regardless. In the case of actual duplicates, the duplicate flag will be ...


6

Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.


6

As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of samtools collate. … | samtools fastq -F 0x900 -@ 48 \ -0 /dev/null -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz - Your samtools fastq command is not doing anything to siphon off ...


6

There are duplicates, in this line: 1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356 According to samtools documentation for flagstat: Provides counts for each of 13 categories based primarily on bit flags in the FLAG field. Each category in the output is broken down into QC pass and QC fail. In the default output format, these are ...


5

You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use. Edit: The FastaVariant class is new in pyfaidx, I hadn't ...


5

Another approach is htsbox. You can get a candidate list with: htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf Here, -q sets min mapping quality, -Q sets min base quality, -v outputs variants only -c outputs VCF, -C gives you base counts on both strands and finally -s5 requires at least 5 high-quality bases to call out an allele. It is ...


5

BioPython has some good tools for processing reads and alignments. http://biopython.org/DIST/docs/tutorial/Tutorial.html There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls https://pysam.readthedocs.io/en/latest/ I would use subprocess to call the aligner and specify the output to a bam ...


5

You might have to manually strip those auxiliary tags off: samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam - The sed bit is searching for the aux tag and removing everything up to the next tab.


5

(1) It won't be super fast but you can provide grep with a file of QNAMES. samtools view file.bam | grep -f 'qnames.txt > subset.sam where qnames.txt has EXAMPLE:QNAME1 EXAMPLE:QNAME2 EXAMPLE:QNAME3 EXAMPLE:QNAME4 EXAMPLE:QNAME5 (2) This would be a little more complicated but can you give an example where the grep might be have the correct QNAME? (3)...


5

Your first command only needs a slight modification to add in -h. This will create a SAM file with a header. samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam -h Include the header in the output. If your primary end goal is to create a BAM, then as Pierre has pointed out you can create BAM files directly from samtools view ...


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

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to: bwa mem (...options) | samtools view -o out.bam samtools fastq -f 4 out.bam > unmatched.fastq


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