NGS one-liner to call variants

These are the variable we later use in the pipeline. We can change [BROKEN LINK: R1] and [BROKEN LINK: R2] and sample for changing the input. We can edit columns to customize the output.

Assembly is used in the [BROKEN LINK: vep] command and if we want to use GRCh37 we need to install the GRCh37 cache and also need to change the reference genome.

BASEDIR=$(dirname "$(realpath "$0")")

THREADS=16
ASSEMBLY="GRCh38"
REFERENCE="$HOME/reference/GRCh38/Homo_sapiens_assembly38.fasta"
TARGET="$HOME/reference/GRCh38/hg38.refGene.exon.bed"
AMPLICON="NO"

R1="$HOME/sample/U0a_CGATGT_L001_R1_005.fastq.gz"
R2="$HOME/sample/U0a_CGATGT_L001_R2_005.fastq.gz"
SAMPLE="U0a"

OUTPUT_DIR="results.$$"

OUTPUT="$OUTPUT_DIR/$SAMPLE"
columns="[%SAMPLE]\t%CHROM\t%POS\t%REF\t%ALT\t%ID\t%FILTER[\t%GT\t%VAF\t%AD\t%DP]\t%Consequence\t%IMPACT\t%SYMBOL\t%Feature\t%EXON\t%INTRON\t%HGVSc\t%HGVSp\t%cDNA_position\t%CDS_position\t%Protein_position\t%Amino_acids\t%Codons\t%Existing_variation\t%MANE_SELECT\t%MANE_PLUS_CLINICAL\t%GENE_PHENO\t%SIFT\t%PolyPhen\t%DOMAINS\t%AF\t%gnomADe_AF\t%gnomADg_AF\t%MAX_AF\t%MAX_AF_POPS\t%CLIN_SIG\t%PHENO\t%PUBMED\t%CANONICAL\n"
header=$(echo "$columns" | sed "s/%//g;s/\[//g;s/\]//g")

mkdir -p "$OUTPUT_DIR"

cp "$0" "$OUTPUT_DIR"

Here we copy the script to result directory to log the command we are running.

Multiqc can show these in a neat table. We just need some formating.

{
│   printf 'oneliner: "%s"\n' "$VERSION";
│   printf 'fastp: "%s"\n' "$(fastp 2>&1 | grep version | cut -d " " -f 2)"
│   printf 'bwa: "%s"\n' "$(bwa 2>&1 | grep Version | cut -d: -f2)"
│   printf 'samtools: "%s"\n' "$(samtools version | sed 1q | cut -d " " -f 2)"
│   printf 'bcftools: "%s"\n' "$(bcftools version | sed 1q | cut -d " " -f 2)"
│   printf 'ensembl-vep: "%s"\n' "$(~/ensembl-vep/vep | grep ensembl-vep | cut -d : -f 2)"
│   printf 'bedtools: "%s"\n' "$(bedtools --version | cut -d " " -f2)"
│   (
│   │   echo "annotation_sources:"
│   │   ~/ensembl-vep/vep --show_cache_info | sed 's/\s/: "/;s/$/"/;s/^/    /'
│   )
} > "$OUTPUT_DIR"/oneliner_mqc_versions.yaml

IGV-reports create easy to browse variant list with alignments. We can add this html to our final html report. At last the multiqc brings all together with custom plots and software versions.

samtools index -@ $THREADS "$OUTPUT.cram"
bcftools index "$OUTPUT.bcf"

. ~/venv/bin/activate

create_report "$OUTPUT".bcf \
│   http://igv-genepattern-org.s3.amazonaws.com/genomes/seq/hg38/hg38.fa \
│   --genome hg38 --flanking 1000 \
│   --sample-columns GT AD DP VAF \
│   --info-columns SYMBOL gnomADg_AF IMPACT Existing_variation \
│   --tracks "$OUTPUT".cram --output "$OUTPUT"_mqc.html
│
samtools stats --reference "$REFERENCE" "$OUTPUT.cram" >"$OUTPUT.cram.stats"
samtools idxstats "$OUTPUT.cram" >"$OUTPUT.cram.idxstats"
samtools flagstat "$OUTPUT.cram" >"$OUTPUT.cram.flagstat"
bcftools stats "$OUTPUT.bcf" >"$OUTPUT.bcf.stats"
"$BASEDIR"/plot_resource_usage.R "$OUTPUT_DIR"
multiqc -f -s -o "$OUTPUT_DIR" "$OUTPUT_DIR"

[BROKEN LINK: pv] is a command line utility that can show a progress bar as data pipe through it. In its manual it uses -s (size) option with `du` command, in our case since our data is compressed we use `gzip -l`. After the first `pv`, it is harder to find out the amount of data passing through so others won’t be accurate but it can still be used to show how much time is passed. Its -S parameter should not be used since it will stop the input prematurely.

They’re commented out because my shell acts weird after the command finishes.

We can start this process in background before the pipeline starts. This gets cpu, memory percentages with `ps` command and file sizes with `du` command every 5 seconds.

monitor_resources() {
│   while ps $$ >/dev/null; do
│   │   du -b $OUTPUT_DIR/* | sed "s#^#$(date +%Y/%m/%d/%H:%M:%S) #" >>"$OUTPUT_DIR/file_sizes.$$.log"
│   │   ps --ppid $$ --forest --no-heading -o %cpu,%mem,cmd 2>/dev/null |
│   │   │   cut -d " " -f 1-6 |
│   │   │   sed "s#^#$(date +%Y/%m/%d/%H:%M:%S) #" |
│   │   │   grep -v "CMD\|pv\|ps" |
│   │   │   awk '{print $4"_"$5"\t"$1"\t"$2"\t"$3}' >>"$OUTPUT_DIR/resources.$$.log"
│   │   sleep 5
│   done
}
monitor_resources &

I would assume no. We have random reads coming from the fastq file and we should wait for alignment process to complete before going on with the variant calling. However, this doesn’t happen because I assume `samtools sort` buffers the whole thing while sorting.

While testing with a smaller subset running in pipe and distinct steps gives the same number of resulting variants. However, when we remove the `samtools sort` I got no variants which in concordance with my initial hunch. This would need more thorough testing to be sure like showing bwa and bcftools running at the same time (but i didn’t keep those stats).

bar@debiantestin:~/runs$ ls -1 results.*/Test.*.bcftools.log | xargs grep '^Lines'

Seperate Steps
results.234225/Test.234225.bcftools.log:Lines   total/split/joined/realigned/skipped:   20/0/0/0/0

Pipe no samtools sort
results.235225/Test.235225.bcftools.log:Lines   total/split/joined/realigned/skipped:   0/0/0/0/0

Pipe
results.70416/Test.70416.bcftools.log:Lines   total/split/joined/realigned/skipped:     20/0/0/0/0

This is not a pipe related concern but after setting up the pipe and running with relatively larger data which is still around 500Mb it took 15 hours.

BWA completes about less than 10 minutes. In this time the first chunk of the alignments are written about 60-70 Mb, then alignments are written incrementally until sort and mpileup finishes.

Investigating further I found we are calling anything and everything; half a million variants in the file. Target intervals only includes 13007 of them which only 7 has depth and quality bigger than 5.

bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -H | wc -l
532845
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -i 'FORMAT/DP>1' -H | wc -l
70674
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -i 'FORMAT/DP>20' -H | wc -l
2408
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -i 'FORMAT/DP>20&&QUAL>20' -H | wc -l
2080
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -i 'FORMAT/DP>20&&QUAL>100' -H | wc -l
1288
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf  -H -T /home/bar/reference/GRCh38/hg38.refGene.exon_padding.bed |wc -l
13007
bar@debiantestin:~/runs/results.89097$ bcftools view U0a.bcf -i 'FORMAT/DP>5&&QUAL>5' -H -T /home/bar/reference/GRCh38/hg38.refGene.exon_padding.bed |wc -l
7

After filtering for target region and the depth and quality we get it down to 15 minutes. For some reason target region overlaps low quality variants.

How much of the variants filtered effects the time it tooks to complete the pipeline. Biggest time concern after eliminating variant count would be the VEP. It takes around ~15 mins to annote 7 variants. Runing the VEP by itself takes only around 15 secs. Meaning there are some other bottlenecks.

We utilize more of the system resources passing around the $THREADS variable but it would be more interesting challenge to scale this pipeline horizontally.

bwa mem -t "$THREADS" -R "@RG\tID:$SAMPLE\tSM:$SAMPLE\tPL:illumina\tLB:lib1\tPU:foo" "$REFERENCE" "$R1" "$R2" 2>"$OUTPUT.$$.bwa.log" > "$OUTPUT.sam"
│   # pv -cN bwa |
samtools collate -@ "$THREADS" -o "$OUTPUT.collated.bam" "$OUTPUT.sam"
samtools fixmate -@ "$THREADS" -m "$OUTPUT.collated.bam" "$OUTPUT.fixmate.bam"
samtools sort -@ "$THREADS" "$OUTPUT.fixmate.bam" -o "$OUTPUT.sorted.bam"  -O bam
samtools markdup -@ "$THREADS" "$OUTPUT.sorted.bam" "$OUTPUT.markdup.bam"
│   # ([ "$AMPLICON" = "NO" ] && samtools markdup -@ "$THREADS" - - || cat) |
│   # pv -cN samtools_markdup |
│
bcftools mpileup --threads "$THREADS" -Ou -A -T "$TARGET" -d 10000 -L 10000 -a "FORMAT/AD,FORMAT/DP" -f "$REFERENCE" "$OUTPUT.markdup.bam" 2>>"$OUTPUT.$$.bcftools.log" |

Even after going through this section I had different number of variants while running the pipeline. It would result in 6, 7 or 8 variants in every ten batches of my test runs. Debuging this was like disassembling an engine and testing each part individually. Variants are counted by [BROKEN LINK: bcftools_norm] step. So whatever was happening before this step. I worked my way up the ladder.

I had 10 different runs with various number of variants. I first make sure [BROKEN LINK: bcftools_mpileup] and [BROKEN LINK: bcftools_call] worked correctly by running it ten times for different number of variant called runs. It produced the same output from every time whether it was 6, 7 or 8 variants. So it was something upstream.

I separated each command before and checked the stats.

for bam in $(find -type f -name "*.*am"); do echo $bam; samtools stats $bam > $bam.stats; done
find -type f -name "*.stats" | xargs grep -Ri "reads mapped and paired"

./results.148485/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.148485/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.148485/U0a.sam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.148485/U0a.collated.bam.stats:SN      reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.148485/U0a.fixmate.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.159805/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.159805/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.159805/U0a.sam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.159805/U0a.collated.bam.stats:SN      reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.159805/U0a.fixmate.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.155247/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229236 # paired-end technology bit set + both mates mapped
./results.155247/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229236 # paired-end technology bit set + both mates mapped
./results.155247/U0a.sam.stats:SN       reads mapped and paired:        4229236 # paired-end technology bit set + both mates mapped
./results.155247/U0a.collated.bam.stats:SN      reads mapped and paired:        4229236 # paired-end technology bit set + both mates mapped
./results.155247/U0a.fixmate.bam.stats:SN       reads mapped and paired:        4229236 # paired-end technology bit set + both mates mapped
./results.150020/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.150020/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.150020/U0a.sam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.150020/U0a.collated.bam.stats:SN      reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.150020/U0a.fixmate.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.162839/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.162839/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.162839/U0a.sam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.162839/U0a.collated.bam.stats:SN      reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.162839/U0a.fixmate.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.156764/U0a.sorted.bam.stats:SN        reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped
./results.156764/U0a.markdup.bam.stats:SN       reads mapped and paired:        4229234 # paired-end technology bit set + both mates mapped

samtools outputs where consistent in itself as it showed same counts in each run but 155247 in particular had totally different counts. Pointing there was something before or during the alignment.

Next I checked the fastq files produced by fastp.

│find -type f -name "*.fq.gz"  | xargs md5sum

a90dd1e71c9e5432e8b67a54e65424a1  ./results.148485/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.159805/U0a.fq.gz
8f716b9f65748efdb3f30ba7d7794d4d  ./results.155247/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.150020/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.162839/U0a.fq.gz
8b4619fd39555cee4d5cca430d999379  ./results.156764/U0a.fq.gz
8b4619fd39555cee4d5cca430d999379  ./results.153722/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.158278/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.151555/U0a.fq.gz
a90dd1e71c9e5432e8b67a54e65424a1  ./results.161314/U0a.fq.gz

It creates three different outputs. Searching the github I found related issues: https://github.com/OpenGene/fastp/issues/506. Removing the fastp step gave the consistent alignment and variant call results.

• Changed the title from NGSoneliner to NGS one-liner to call variants and made additions in order to make the subject clear thanks to GenoMax’s suggestions.