|
| 1 | +# 2. Quality control {.unnumbered} |
| 2 | + |
| 3 | +::: callout-warning |
| 4 | +# Important! |
| 5 | + |
| 6 | +In the next steps we are going to copy-paste code, adjust it to our needs, and execute it on the command-line. |
| 7 | + |
| 8 | +**Please open a plain text editor to paste the code from the next steps, to keep track of your progress!** |
| 9 | +::: |
| 10 | + |
| 11 | +For simplicity's sake, most steps will be geared towards an analysis of a single sample. It is recommended to follow a basic file structure like the following below: |
| 12 | + |
| 13 | +``` |
| 14 | +my_project/ |
| 15 | +├── raw_data/ # Contains barcode directories |
| 16 | +│ ├── barcode01 # Contains the raw, gzipped FASTQ files |
| 17 | +│ ├── file1.fastq.gz |
| 18 | +│ └── file2.fastq.gz |
| 19 | +│ ├── barcode02 |
| 20 | +│ └── file1.fastq.gz |
| 21 | +│ └── barcode03 |
| 22 | +│ ├── file1.fastq.gz |
| 23 | +│ ├── file2.fastq.gz |
| 24 | +│ └── file3.fastq.gz |
| 25 | +└── results/ # This is where the output files will be stored |
| 26 | +└── log/ # This is where log files will be stored |
| 27 | +``` |
| 28 | + |
| 29 | +When running any command that generates output files, it's essential to ensure that the output directory exists *before* executing the command. While some tools will automatically create the output directory if it's not present, this behavior is not guaranteed. If the output directory doesn't exist and the tool doesn't create it, the command will likely fail with an error message (or, worse, it might fail silently, leading to unexpected results). This is not required if you are running a snakemake workflow. |
| 30 | + |
| 31 | +To prevent a lot of future frustration, create your output directories beforehand with the `mkdir` command as such: |
| 32 | + |
| 33 | +``` bash |
| 34 | +mkdir -p results |
| 35 | +mkdir -p log |
| 36 | +# Create a subdirectory |
| 37 | +mkdir -p results/assembly |
| 38 | +``` |
| 39 | + |
| 40 | +To use the required tools, activate the Singularity container as follows: |
| 41 | +```bash |
| 42 | +singularity shell naam_workflow.sif |
| 43 | +``` |
| 44 | + |
| 45 | +## 2.1 Merging and decompressing FASTQ {.unnumbered} |
| 46 | + |
| 47 | +Any file in linux can be pasted to another file using the cat command. zcat in addition also unzips gzipped files (e.g. .fastq.gz extension). If your files are already unzipped, use cat instead. |
| 48 | + |
| 49 | +**Modify and run:** |
| 50 | + |
| 51 | +``` bash |
| 52 | +zcat {input.folder}/*.fastq.gz > {output} |
| 53 | +``` |
| 54 | + |
| 55 | +- `{input.folder}` should contain all your .fastq.gz files for a single barcode. |
| 56 | +- `{output}` should be the name of the combined unzipped fastq file (e.g. all_barcode01.fastq). |
| 57 | + |
| 58 | +## 2.2 Running fastp quality control software {.unnumbered} |
| 59 | + |
| 60 | +The [fastp](https://github.com/OpenGene/fastp){target="_blank"} software is a very fast multipurpose quality control software to perform quality and sequence adapter trimming for Illumina short-read and Nanopore long-read data. |
| 61 | + |
| 62 | +Because we are processing Nanopore data, several quality control options have to be disabled. The only requirement we set is a minimum median phred quality score of the read of 10 and a minimum length of around the size of the amplicon (e.g. 400 nucleotides). |
| 63 | + |
| 64 | +```bash |
| 65 | +fastp -i {input} -o {output} -j /dev/null -h {report} \ |
| 66 | +--disable_trim_poly_g \ |
| 67 | +--disable_adapter_trimming \ |
| 68 | +--qualified_quality_phred 10 \ |
| 69 | +--unqualified_percent_limit 50 \ |
| 70 | +--length_required {min_length} \ |
| 71 | +-w {threads} |
| 72 | +``` |
| 73 | + |
| 74 | +- `{input}` is the merged file from step 2.1. |
| 75 | +- `{output}` is the the quality controlled `.fastq` filename (e.g. `all_barcode01_QC.fastq`). |
| 76 | +- `{report}` is the QC report filename, containing various details about the quality of the data before and after processing. |
| 77 | +- `{min_length}` is the expected size of your amplicons, to remove very short "rubbish" reads, generally the advise is to set it a bit lower than the expected size. Based on the QC report, which lists the number of removed reads you may adjust this setting, if too many reads are removed. |
| 78 | + |
| 79 | +::: callout-note |
| 80 | +`{threads}` is a recurring setting for the number of CPUs to use for the processing. On a laptop this will be less (e.g. 8), on an HPC you may be able to use 64 or more CPUs for processing. However, how much performance increase you get depends on the software. |
| 81 | +::: |
| 82 | + |
| 83 | +## 2.3 Mapping reads to primer reference {.unnumbered} |
| 84 | + |
| 85 | +To precisely trim the primers we map the reads to a reference sequence based on which the primers were designed. This is to make sure, when looking for the primer locations, all primer location can be found. To map the reads we use [minimap2](https://github.com/lh3/minimap2) with the `-x map-ont` option for ONT reads. `-Y` ensures reads are not hardclipped. Afterwards we use [samtools](https://www.htslib.org/) to reduce the `.bam` (mapping) file to only those reads that mapped to the reference and sort the reads in mapping file based on mapping position, which is necessary to continue working with the file. |
| 86 | + |
| 87 | +```bash |
| 88 | +minimap2 -Y -t {threads} -x map-ont -a {reference} {input} | \ |
| 89 | +samtools view -bF 4 - | samtools sort -@ {threads} - > {output} |
| 90 | +``` |
| 91 | + |
| 92 | +- `{reference}` is the fasta file containing the reference that your primers should be able to map to. |
| 93 | +- `{input}` is the QC fastq file from step 2.2. |
| 94 | +- `{output}` is the mapping file, it could be named something like `barcode01_QCmapped.bam` |
| 95 | + |
| 96 | +## 2.4 Trimming primers using Ampliclip {.unnumbered} |
| 97 | + |
| 98 | +[Ampliclip](https://github.com/dnieuw/Ampliclip) is a tool written by [David](https://github.com/dnieuw/) to remove the primer sequences of nanopore amplicon reads. It works by mapping the primer sequences to a reference genome to find their location. Then it clips the reads mapped to the same reference (which we did in the previous step) by finding overlap between the primer location and the read ends. It allows for some "junk" in front of the primer location with `--padding` and mismatches between primer and reference `--mismatch`. After clipping it trims the reads and outputs a clipped `.bam` file and a trimmed `.fastq` file. `--minlength` can be set to remove any reads that, after trimming, have become shorter than this length. Set this to the value that was used in the QC section (e.g. 400). |
| 99 | + |
| 100 | +After the trimming the clipped mapping file has to be sorted again. |
| 101 | + |
| 102 | +```bash |
| 103 | +samtools index {input.mapped} |
| 104 | + |
| 105 | +ampliclip \ |
| 106 | +--infile {input.mapped} \ |
| 107 | +--outfile {output.clipped}_ \ |
| 108 | +--outfastq {output.trimmed} \ |
| 109 | +--primerfile {primers} \ |
| 110 | +--referencefile {reference}\ |
| 111 | +-fwd LEFT -rev RIGHT \ |
| 112 | +--padding 20 --mismatch 2 --minlength {min_length} > {log} 2>&1 |
| 113 | + |
| 114 | +samtools sort {output.clipped}_ > {output.clipped} |
| 115 | +rm {output.clipped}_ |
| 116 | +``` |
| 117 | + |
| 118 | +- `{input.mapped}` is the mapping file created in step 2.3. |
| 119 | +- `{output.clipped}` is the mapping file processed to clip the primer sequences off (e.g. `barcode01_clipped.bam`). |
| 120 | +- `{output.trimmed}` is the trimmed fastq file, this contains all reads mapped to the reference with primer sequences trimmed off (e.g. `barcode01_trimmed.bam`). |
| 121 | +- `{primers}` is the name of the primer sequence fasta file. Make sure names of the primers have either 'LEFT' or 'RIGHT' in their name to specify if it is a left or right side primer. |
| 122 | +- `{reference}` is the name of the reference file, this must be the same file as was used for mapping in section 2.3. |
| 123 | +- `{min_length}` is the minimum required length of the trimmed reads, set it to the same value as when using `fastp`. |
| 124 | + |
| 125 | +To see what has happened in the trimming process we can open the `.bam` mapping files before and after primer trimming using the visualization tool [UGENE](https://ugene.net/), a free and open source version of the software [geneious](https://www.geneious.com/). |
| 126 | + |
| 127 | +In UGENE you can open a `.bam` via the "open file" option. |
| 128 | + |
| 129 | +::: {.callout-note} |
| 130 | +We now have our quality controlled sequence reads which we can use to create a consensus sequence in the next chapter. |
| 131 | +::: |
0 commit comments