Hot answers tagged

9

samtools quickcheck is all you need. From the manual: Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only). Data in the middle of the file ...


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


4

Usually it doesn't actually matter, but if you want it to be very correct: RGID=SomeLibraryID \ RGLB=SomeLibraryID \ RGPL=illumina \ RGPU=Some_Sequencer_ID \ RGSM=SomeSampleID That assumes you have one library per sample and aren't splitting the output BAM files by lane.


2

GATK cannot combine VCF files generated by Sniffles, the variant caller that I used to call structural variants (answer from the GATK team).


2

Answered in the GATK forum Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.


2

It sounds like you want samtools faidx foo.fa followed by samtools faidx foo.fa chr1 > your_subset_file.fa (or whatever the first chromosome is). The output file is then a regular fasta file subset as you like. You can get any chromosome you want with that. In fact, you can also do regions (e.g., chr1:100-1000, though note that the sequence name in the ...


2

Alternatively, I really like using bbduk which is part of the BBMap suite. I've processed every nascent sequencing dataset that has been published, and found a lot of quirky errors with older datasets using TrimGalore. bbduk is a little more fine-tuneable relative to cutadapt/trimmomatic/trimGalore (built on top of cutadapt)/fastp and the run time and ...


2

You're best off just using fastp or Trim Galore!, both of which will determine the adapter sequence for you. Trim Galore! uses a built-in list of known sequences for this, whereas fastp uses read overlap.


2

@StupidWolf's answer is correct -- that first number in the flagstat output is what you want to look at to see the number of reads marked as duplicates. I wanted to add that the number given in the Picard metrics file is in fact 6%, not 0.06%. This is due to an ambiguity that is widespread in Picard metrics files; in many places the program emits a fraction ...


2

From the Picard documentation: DUPLICATION METRICS: Metrics that are calculated during the process of marking duplicates within a stream of SAMRecords. UNMAPPED_READS The total number of unmapped reads examined. (Primary, non-supplemental) It won't alter the flags on these reads, but it will count them in the summary report it generates. You should be ...


1

I'm late to the party but simply run cellranger mkfastq twice, with different arguments for the mask and with --filter-dual-index for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace ...


1

Picard, like most tools, prints logs to stderr. It should be in the Galaxy history. Are you using a version of picard from the toolshed (if so, which one)? You might set cleanup_job to never in the config, which will prevent galaxy from cleaning the working directory. This will enable seeing what the job actually produced. Most likely there's more in the ...


1

As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.


1

The proper solution is to use a different tool (some of this you could do with computeGCBias in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s): $ cat foo.awk BEGIN{split(excludes, excludeList, " ")} {exclude=0 if($1 == "@SQ") { for(ex in excludeList) { if(ex == substr($2, ...


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