In this report you will find all the information necessary to follow the steps to analyze SARS-CoV-2 data with Galaxy.
During this training we will following these steps:
- Register/Login: Register or login into Galaxy website
- Data: Upload data for the analysis
- Quality: Analysis of the quality of the raw reads
- Trimming: Quality trimming using fastp
- Mapping: Mapping reads to reference genome with Bowtie2
- Stats: Mapping statistics with picard
- Amplicons: Preprocessing steps mandatory for amplicon sequencing data
- Variants: Variant calling and annotation
- Consensus: Consensus genome generation
- Additional info: Websites for additional info.
First of all, you have to create a user (or login into your user if you already have one) in the European Galaxy website by clicking in the button indicated in the image:
You can do this analysis without creating a user, but your history will be saved into your user so you can access later at home.
Once you are registered we can start with the analysis workflow.
From now on, each job we run in Galaxy will have a unique number for identifying each process. This numbers can differ depending on the number of samples and the times you run or delete any process. This training's snapshots were taken using two samples and some process were deleted for any reason, so numbers MAY DIFFER. Also, this exercise's snapshots were taken with other dataset, so sample names may differ too.
Before starting with any analysis we have to upload de data we want to analyze into galaxy. The are three different ways to upload data to Galaxy, from which we will explain you some of them along this training.
In order to do upload files from URL we have to follow these steps:
- In the left side panel, select Upload
- In the new panel select Paste/Fetch Data
- Then copy the following block of text for the samples:
https://zenodo.org/record/5724464/files/SARSCOV2-1_R1.fastq.gz?download=1
https://zenodo.org/record/5724464/files/SARSCOV2-1_R2.fastq.gz?download=1- Now, in the Download data from the web by entering URLs (one per line) or directly paste content. square, paste the text you copied before
- Select Start
- When everything is green in the screen, select Close
Once files are uploaded we can close the window clicking on "close" button. Both files will appear in the right panel.
If we select the "eye" in the datasets from the history, we can see the content of the data:
Once we have the raw data, an important step is to analyze the quality of the reads, to know if the reads are worth it. To do this, we have to:
- Search for the fastqc tool
- Select FastQC Read Quality reports and set the following parameters:
- Select Multiple datasets in Raw read data from your current history
- then select in the bar the two datasets from the history
- Then go down and select Run tool
This program will generate a message like this one, were we can read that, each .fastq file is going to generate two different jobs, one for the Raw Data and another one for the .HTML report.
To visualize the information coming from FastQC we just have to select the job of interest. In this case we are interested in the "Web page results" so for the sample we want to see the results we have to click in the eye to visualize galaxy results:
This report gives different type of information about the sequencing quality of the reads. Clicking in the arrows in the bottom right/left corners of the page we can hide/show the side panels.
So the central panel with the results we want to visualize will bee better seen.
For more information about FastQC output visit FastQC website
Once we have check the quality of our reads, it's important to trim low quality nucleotides from those reads, for which we will use Fastp. So, in the search bar you look for fastp and then select "fastp - fast all-in-one preprocessing for FASTQ files". There, we will have to change some parameters ensure the trimming accuracy for this amplicon data. First of all we are going to do the analysis for the sample we gave to you (201569). These are the field we will have to change:
3.Single-end or paired reads > Paired
4.Input 1 > Select SARSCOV2-1_R1.fastq.gz?download=1
Input 2 > Select SARSCOV2-1_R2.fastq.gz?download=1
5.Display Filter Options
-
Quality Filtering options
6.Qualified Quality Phred = 20
7.Unqualified percent limit = 10
-
Length Filtering Options
8.Length required = 50
9.Read modification options
10.PoliX tail trimming > Enable polyX tail trimming
-
Per read cutting by quality options
11.Cut by quality in front (5') > Yes
12.Cut by quality in tail (3') > Yes
13."Execute"
A message like this one will appear, which means that 4 results will be generated:
- One with the R1 trimmed reads:
fastp on data 2 and data 1: Read 1 output - Another one with the R2 trimmed reads:
fastp on data 2 and data 1: Read 2 output - Another one with the HTML results:
fastp on data 2 and data 1: HTML report - A last one which is not relevant for this training
Once fastp analysis is done, you can see the results by clicking in the eye ("View Data") in the fatp HTML results. We will see a report like this one:
This report can be downloaded by clicking in the results name ("fastp on data 2 and data 1: HTML report") and then in "Download".
This will download a .zip folder that contains the .html report, that can be visualize in your computer using any browser (such as Google Chrome) and interactively scroll around the results.
Among the most relevant results, you have the:
- Summary: Stats summary
- After filtering: Statistics of the reads after quality filtering
- reads passed filters: Reads remaining after quality filter trimming
- reads with low quality: Reads that were remove due to low quality
- reads too short: Reads that didn't pass the minimum length filter.
- After filtering: Plots after filtering
- After filtering: read1: quality: Plot with the evolution of R1 quality over read position. Usually it decays in the last nucleotides.
- After filtering: read2: quality: Same plot for R2.
- After filtering: Statistics of the reads after quality filtering
For more information about FastQC output visit Fastp github
In order to call for variants between the samples and the reference, it's mandatory to map the sample reads to the reference genome. To do this we need the fasta file of the reference and the Bowtie2 index of that fasta file.
Now we can start with the main mapping process. The first thing we have to do is look for the program "Bowtie2" in the search bar and then select "Bowtie2 - map reads against reference genome". Here we will have to set the following parameters, for the first sample:
3.Is this single or paired library > Paired-end
4.Fasta/Q file #1: fastp on data 2 and data 1: Read 1 output
5.Fasta/Q file #2: fastp on data 2 and data 1: Read 2 output
-Will you select a reference genome from your history or use a built-in index?
6.Select reference genome > `SARS-CoV-2 isolate Wuhan-Hu-1, complete genome (NC_045512.2)`
7.Do you want to use presets? > Very sensitive local
8.Save the bowtie2 mapping statistics to the history > Yes
9.Run tool
We will see a message like this one:
Now we can see the mapping results for the samples. The bowtie2 resulting file is a .bam file, which is not easy to read by humans. This .bam file can be downloaded by clicking in the alignment file and then into download. Then, the resulting .gz file will contain the alignment .bam file that can be introduced in a software such as IGV with the reference genome fasta file.
In our case, the file that can be visualize is the statistics file, which contains information such as the percentage of reads that aligned.
The previously shown files give few human readable information, because mapping files are supposed to be used by other programs. In this sense, we can use some programs to extract relevant statistical information about the mapping process.
Another program that gives statistical information about the mapping process is Picard. To run this program you just have to search "Collect Wgs Metrics" and then select "CollectWgsMetrics compute metrics for evaluating of whole genome sequencing experiments".
- Using reference genome >
SARS-CoV-2 isolate Wuhan-Hu-1, complete genome (NC_045512.2) - Run tool
This process will generate one output file per .bam alignment file selected as input.
Picard results consist in quite long files, so the best is to download those results and visualize them in your computer. Yo you have to click in the CollectWgsMetrics job you want to download, and then click in the save button:
Then you just have to open the file with Excell in your computer, and you will see a file with different columns with information about the percentage of the reference genome that is covered by the reads at a specific depth or the mean depth of coverage of the reference genome.
So in the results table you can see that the "Mean Coverage" is around 115, which means that each position in the reference genome is supported by 115 reads, by mean. The "PCT 10X" represents the percentage of the reference genome that is sequence with a 10X depth, which represents the percentage of the reference genome that is at least sequenced 10 times (79,57%).
After mapping the reads to the reference genome, we are interested in removing the sequences of the amplicon primers. To do that you will use a program called iVar, and you will need a bed file with the positions of those amplicon primers.
You have to search for "ivar trim" in the search bar and select "ivar trim Trim reads in aligned BAM". Then follow these steps:
- Source of primer information >
Built-in - Primer scheme name >
SARS-CoV-2 ARTICv3 - Include reads not ending in any primer binding sites? >
Yes. - Require a minimum length for reads to retain them after any trimming? >
Yes, and provide a custom threshold - Minimum trimmed length threshold =
20 - Run tool
The resulting file from iVar will be a new BAM file where amplicon primer positions will be removed, so there's no result to visualize.
Once we have the alingment statistics and files with amplicon primers trimmed, we can start with the variant calling process.
We have to upload the fasta reference genome, and we are going to do it using the following URL:
https://zenodo.org/record/5724970/files/GCF_009858895.2_ASM985889v3_genomic.200409.fna.gz?download=1So you have to copy the URL direction, and do then select "Download from web or upload from disk", the select "Paste/Fetch data" and in the new window that appears you should paste the URL above and finally select "Start". Now the fasta file of the reference is download.
To call for variants between the sample and the reference we are going to use iVar variants, so you have to search for "ivar" in the search bar and select "ivar variants Call variants from aligned BAM file". Then select the following parameters:
- In "Bam file" make sure it is selected ivar Trimmed bam files
- Reference >
GCF_009858895.2_ASM985889v3_genomic.200409.fna.gz?download=1 - Minimum frequency threshold >
0.25 - Output format >
Both Tabular and VCF - Run tool
iVar will create two output files per each bam file uploaded, one with the tabular results, and another one with the vcf files.
To display iVar variants results, select the 👁️ icon in the right pannel of the Tabular results and then extend the central panel by clicking in the arrows in the bottom.
iVar results consist in a Tab separated file containing all the variants found between the reference and the sample with an Alle Frequency higher than threshold (0.25). Each line represents a variant and the columns give information about that variant, such as the position in the reference genome, the reference allele, the alternate allele, if that variant passed the filters, and so on.
This variants have passed a filter for the minimum quality of the variant, which we set as 20, and the allele frequency threshold. However, we will have a closer look only to those variants present in an allele frequency higher than 0.75 which are the ones that are going to be included in the consensus.
Once we have the variants called, it's interesting to annotate those variants, for which you will use SnpEff. Search for "snpeff" in the searh bar and select "SnpEff eff: annotate variants for SARS-CoV-2", then change the following parameters:
- In "Sequence changes (SNPs, MNPs, InDels)" make sure it is selected iVar variants results
- Run tool
The SnpEff gives three different results, from which the most interesting ones are:
- Snpeff eff: Which is a VCF file with the annotation results. It is a very similar file to the ones we saw before for VarScan and Bcftools but with the last column different, containing relevant information about that variant.
- Snpeff eff HTML stats: This file is a CSV file that contains statistics about the variant annotation process, such as the percentage of variants annotated, the percentage of variants that are MISSENSE or SILENT, the percentage that have HIGH, MODERATE or LOW effect, and so on.
Once we have the most relevant variants that can be considered to include in the consensus genome, you can start with the consensus genome generation.
Now we are going to generate the consensus genome using iVar. We are going to search for ivar and the select "ivar consensus Call consensus from aligned BAM file". Now follow these steps:
- In "BAM fle" make sure iVar trim bam is selected
- Minimum frequency threshold >
0.75 - Run tool
This will just generate a fasta file identical to the reference one, except for those nucleotides that are variants from the VCF file.
Now we are going to determine the lineage of the samples. We will use a software called pangolin. We are going to use the masked consensus genomes generated in the previous steps as follows:
- Search for the pangolin tool
- Select Pangolin Phylogenetic Assignment of Outbreak Lineages and set the following parameters:
- In "Input FASTA File(s)" make sure iVar consensus fasta is selected
- Version of pangolin-data to use >
Download latest pangolin-data version from web - Version of constellations to use >
Download latest pangolin-data version from web - Run tool
Now we are going to have a look to the results from pangolin:
As you can see, results are in table format, where you have in first place the reference genome, then de lineage assingned, which seems to be B.1.1.7 lineage. This lineage corresponds as you can see in the 5th column to the Alpha lineage.
If you have any problem following this training and you want to visualize the resulting file you can access them through this URL
And viralrecon workflfow in this URL
- GiSaid: The GISAID Initiative promotes the rapid sharing of data from all influenza viruses and the coronavirus causing COVID-19. This includes genetic sequence and related clinical and epidemiological data associated with human viruses, and geographical as well as species-specific data associated with avian and other animal viruses, to help researchers understand how viruses evolve and spread during epidemics and pandemics.
- ENA: The European Nucleotide Archive (ENA) provides a comprehensive record of the world’s nucleotide sequencing information, covering raw sequencing data, sequence assembly information and functional annotation.
- Nextstrain: Nextstrain is a project to harness the scientific and public health potential of pathogen genome data. Our goal is to aid epidemiological understanding of pathogen spread and evolution and improve outbreak response.
- outbreak.info: Outbreak.info aggregates data across scientific sources, providing tools to meet three major aims: (1) Track daily developments in SARS-CoV-2 variants. (2) Integrate publications, preprints, clinical trials, datasets, protocols, and other resources into one searchable library of COVID-19 research. (3) Track trends in COVID-19 cases and deaths.