6.1. Effect of the trimming on SNPs detection

Description:Effect of trimming (or not trimming) on the the SNPs detection.

In this case example, we will take a paired-end data set, and apply the quality pipeline using trimming quality (removing bases with quality below 30). Then, we will run the variant calling pipeline to perform the mapping on a reference and detect SNPs.

We will repeat this analysis without trimming low quality reads at all.

We will finally compare the two sets of SNPs showing that the trimming quality is not important in this example. Meaning that the mapping tool used (freebayes) is able to cope with low quality reads.

6.1.1. The data

We will use a paired-end data set (MiSeq 250bp). It contains 250,000 reads (X2). The organism sequenced is Bordetella. As a reference, we use the ENA accession CP010347.1. The data will be posted later but the original data were generated at Pole Biomics (Institut Pasteur) and named Tohama-R0_S4_L001_R1_001 from which we used only the first 250,000 reads.

Here is a boxplot of the base quality across the reads showing that the quality is quite high and falls below 30 after 200 bases.


6.1.2. Quality pipeline

Assuming DATA (fastq.gz files) are in <DIR1> directory, type this command to create the quality pipeline and config file automatically:

sequana --pipeline quality --input-dir <DIR1> --project trimming

Then go to the project and execute the pipeline:

cd trimming
snakemake -s quality.rules --stats report/stats.txt -p -j 4 --forceall

The final cleaned reads are in trimming/report/ (refered to <DIR2> hereafter) and named after the project: (trimming_R1_.cutadapt.fastq.gz and trimming_R2_.cutadapt.fastq.gz). These two files should be used later as the input of the variant_calling pipeline, as shown hereafter.

There is no adapters in the data so in the config file, the adapter sections are empty (no forward or reverse adapters). Note, however, that bad quality bases below 30 (default) are removed. In order to set the quality to another values, use sequana with the –quality option

See also

See the Tutorial and User guide and reference sections for more details.

6.1.3. Variant analysis

The output of the quality pipeline will be the input of the variant calling pipeline:

sequana --pipeline variant_calling --input-dir <DIR2> --project variant_trimming

Here you need to make sure that the config.yaml configuration file has the correct reference. See the Tutorial section (variant section).

reference = "CP010347"
from bioservices import EUtils
eu = EUtils()
data = eu.EFetch(db="nuccore",id=reference, rettype="gbwithparts",
with open("data.gbk", "w") as fout:
from bioservices import ENA
ena = ENA()
data = ena.get_data(reference', 'fasta')
with open("data.fa", "w") as fout:
from sequana import snpeff
v = snpeff.SnpEff("data.gbk")

Edit the config.yaml to change those sections:

# snpEff parameter
    do: yes
    reference: "data.gbk"

# Bwa parameter for reference mapping
    reference: "data.fa"

Run the analysis:

cd variant_trimming
snakemake -s variant_calling.rules --stats report/stats.txt -p -j 4 --forceall

Once done, you should have VCF files in variant/report/ named cutadapt.ann.vcf

6.1.4. No trimming

Repeat the previous two steps. In the first step, change the adapter section (cutadapt) to set the quality to zero (this prevents the trimming of bad quality bases):

    quality: 0,0

Change the project name e.g. no_trimming as a tag to the project in the first step and variant_no_trimming.

6.1.5. SNPs results comparison

You should now have two VCF files. Here below we plot the read depth versus strand balance. The color will indicates the overall freebayes score (normalised by the largest score). A good candidate should have large score and balance value around 0.5 The y-axis shows the read depth.

from pylab import *

from sequana import vcf_filter

vcf1 = vcf_filter.VCF("variant/report/cutadapt.ann.vcf")
vcf2 = vcf_filter.VCF("variant_no_trimming/report/variant_no_trimming.ann.vcf")
    df1 = vcf1.vcf_to_csv("dummy")
df2 = vcf2.vcf_to_csv("dummy")

scatter(list(df1.strand_balance.values), list(df1.depth.values),
xlabel("strand balance")

scatter(list(df2.strand_balance.values), list(df2.depth.values),
title("Trimming quality (left) vs no trimming (right)

In this figure the LHS (trimming) 294 SNPs were found while in the RHS (no trimming) 309 were found. The additional SNPs all have low coverage below 20. A third of them have low balance strand.

There is one SNP found in the trim case not found in no_trim. However, it is marginal with strand balance of 0.12, depth of 11, frequence of 0.73 and one of the lowest score

6.1.6. Conclusions

The detection of SNPs does not suffer from not trimming low quality bases below 30. Actually, some new SNPs are found. However, the are usually not significant (low depth, low score or unbalanced). Interestingly, the distribution of the SNPs in the depth vs strand balance plane seems to be more centered on strand balance=0.5 could be interesting to extend the analysis to more data, lower quality, or higher quality threshold.