4

Going VCF to mpileup is not really something one does or can do. The mpileup should be generated from a BAM or SAM file or something else that has raw, unfiltered read alignments in it. The VCF just has data on variant sites, and usually just variant sites that passed a certain likelihood threshold at that. With the VCF, the best you can do is generate a ...


3

Instead of doing this kind of tweaking I would write my script in a way that these options are passed as strings that can be empty. Something like: #/ dummy syntax just for illustration: if(you_want_filtering){ filter_f='-f something' filter_F='-F something_else' } else { filter_f=filter_F='' } samtools view (...) ${filter_f} ${filter_F} So if you ...


3

It seems you do not have bowtie installed or loaded. Use module load bowtie2 together when you load samtools.


2

Use samtools flagstat with option -O tsv: Using -O tsv selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names. SEE ALSO: What Does Samtools Flagstat ...


2

Htslib is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM: #include "htslib/sam.h" int main(int argc, char *argv[]) { samFile *fp = sam_open(argv[1], "r"); bam_hdr_t *h = sam_hdr_read(fp); bam1_t *b = bam_init1(); while (sam_read1(fp, h, b) >= 0) { if (...


2

#include <stdio.h> #include "htslib/sam.h" int main(int argc, char *argv[]) { samFile *fp = sam_open(argv[1], "r"); hts_idx_t *idx = sam_index_load(fp, argv[1]); bam_hdr_t *h = sam_hdr_read(fp); bam1_t *b = bam_init1(); hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]); while (sam_itr_next(fp, ...


2

Using GNU parallel: ls *.bam | parallel -j 2 "samtools fastq {} > {.}.fastq" Use -j to control the number of parallel jobs. That only works for single-end data though as paired-end data must first be sorted or collated by name.


2

I would start with: module avail Then look to see (a) if there is a samtools module avilable and (b) exactly what it's called. If it's called samtools/1.11 then: module load samtools/1.11 Should work, if it doesn't speak to your sys admin.


1

What Devon said. Alternatively, you could also just convert the bam to sam and then look for lines where the 3rd field (the reference the sequence is aligned to) is what you want: samtools view pseudoalignments.bam | awk '$3=="ENST00000367969.8"' But really, it's far better to sort and index the file as Devon suggested so you can use dedicated ...


1

Sort the BAM file again with samtools sort and then run samtools index on the result. That should fix the issue.


1

Related question: Access base aligned to particular reference position After an initial effort posted here, I ended up rewriting this and testing it a little, and posted the script here. It still has some drawbacks (handling indels intelligently) but now it more or less works.


1

If you merge a lot of BAM files you lose the overhead of the header which depending on the size of your BAM files this can be a significant difference or not. With FASTQ files there should be less difference (as they don't have an overall header). However, the compression ratio can change depending on how you zip your files. [EDIT: I just tested it merging ...


1

The helpdesk for my institution's supercomputers showed me the following flags: g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.


1

I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the SeqAn BAM tutorial, it looks like you need to use the BamFileIn object for reading BAM/SAM files: #include <seqan/bam_io.h> using namespace seqan; int main(){ CharString bamFileName = "copyPolly_SRR6511930.bam"; // Open ...


1

If it's just a comparison of two sequences that you're looking for, then diffseq from the emboss toolkit should be okay: https://www.bioinformatics.nl/cgi-bin/emboss/diffseq I think this tool requires the sequences to be oriented in the same direction. If you discover a lot more variants than you expect, try reverse-complementing one of the sequences.


1

My best guess would be a corruption during download as: samtools quickcheck -qvvvv ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam verbosity set to 4 checking ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/...


1

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021. What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the ...


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