RosaBioinfo

RNASeq pipeline of alignment of Culex ovaries metatranscriptomic reads to wPipPel and pWCP plasmid with AliNe

Maxime Mahout

2025-02-18

RNASeq pipeline for alignment of reads
UMI-tools deduplication
Final samtools treatment

In this notebook we want to perform alignment of our raw reads from Culex ovaries metatranscriptomes on the reference genome of the plasmid pWCP.

RNASeq pipeline for alignment of reads

AliNe is a workflow method by Jacques Dainat for performing read trimming, quality assessment and alignment over a reference, useful in particular for RNA-Seq mapping.

We want to utilize 30 high quality metatranscriptomic pairs of reads of Culex quinquefasciatus ovaries from the lab strain MTQ NG (after three days of emergence, and collected at various days after: J1 (Day 1), J3 (Day 3), J7 (Day 7) and J14 (Day 14)).

Our reference wPipPel_Plasmid.fasta is composed of the two reference genomes wPipPel (/Ref_genome/wPipPel_index/GCF_000073005.1/GCF_000073005.1_ASM7300v1_genomic.fna) and the plasmid pWCP (pWCP.fasta, MAG O07 from Reveillaud et al, 2019).

This is important as wPipPel contains copies of IS110 transposase and homologues of toxins so mapping to both allows us to resolve potential ambiguity and non-specific read mapping.

We separate our workflow into a data folder for the data and a references folder for the references.

>ls

data    references

>ls data

## Cxq_NG_J01_ov_10_R1.fastq.gz         Cxq_NG_J01_ov_10_R2.fastq.gz
## Cxq_NG_J01_ov_11_R1.fastq.gz         Cxq_NG_J01_ov_11_R2.fastq.gz
## Cxq_NG_J01_ov_13_R1.fastq.gz         Cxq_NG_J01_ov_13_R2.fastq.gz
## Cxq_NG_J01_ov_14_R1.fastq.gz         Cxq_NG_J01_ov_14_R2.fastq.gz
## Cxq_NG_J03_ov_01_R1.fastq.gz         Cxq_NG_J03_ov_01_R2.fastq.gz
## Cxq_NG_J03_ov_02_R1.fastq.gz         Cxq_NG_J03_ov_02_R2.fastq.gz
## Cxq_NG_J03_ov_03_R1.fastq.gz         Cxq_NG_J03_ov_03_R2.fastq.gz
## Cxq_NG_J03_ov_04_R1.fastq.gz         Cxq_NG_J03_ov_04_R2.fastq.gz
## Cxq_NG_J03_ov_06_R1.fastq.gz         Cxq_NG_J03_ov_06_R2.fastq.gz
## Cxq_NG_J03_ov_07_R1.fastq.gz         Cxq_NG_J03_ov_07_R2.fastq.gz
## Cxq_NG_J03_ov_11_R1.fastq.gz         Cxq_NG_J03_ov_11_R2.fastq.gz
## Cxq_NG_J03_ov_12_R1.fastq.gz         Cxq_NG_J03_ov_12_R2.fastq.gz
## Cxq_NG_J07ns_ov_01_R1.fastq.gz           Cxq_NG_J07ns_ov_01_R2.fastq.gz
## Cxq_NG_J07ns_ov_02_R1.fastq.gz           Cxq_NG_J07ns_ov_02_R2.fastq.gz
## Cxq_NG_J07ns_ov_03_R1.fastq.gz           Cxq_NG_J07ns_ov_03_R2.fastq.gz
## Cxq_NG_J07ns_ov_04_R1.fastq.gz           Cxq_NG_J07ns_ov_04_R2.fastq.gz
## Cxq_NG_J07ns_ov_06_R1.fastq.gz           Cxq_NG_J07ns_ov_06_R2.fastq.gz
## Cxq_NG_J07ns_ov_07_R1.fastq.gz           Cxq_NG_J07ns_ov_07_R2.fastq.gz
## Cxq_NG_J07ns_ov_08_R1.fastq.gz           Cxq_NG_J07ns_ov_08_R2.fastq.gz
## Cxq_NG_J07ns_ov_09_R1.fastq.gz           Cxq_NG_J07ns_ov_09_R2.fastq.gz
## Cxq_NG_J07s_ov_01_R1.fastq.gz            Cxq_NG_J07s_ov_01_R2.fastq.gz
## Cxq_NG_J07s_ov_02_R1.fastq.gz            Cxq_NG_J07s_ov_02_R2.fastq.gz
## Cxq_NG_J07s_ov_03_R1.fastq.gz            Cxq_NG_J07s_ov_03_R2.fastq.gz
## Cxq_NG_J07s_ov_04_R1.fastq.gz            Cxq_NG_J07s_ov_04_R2.fastq.gz
## Cxq_NG_J07s_ov_05_R1.fastq.gz            Cxq_NG_J07s_ov_05_R2.fastq.gz
## Cxq_NG_J07s_ov_07_R1.fastq.gz            Cxq_NG_J07s_ov_07_R2.fastq.gz
## Cxq_NG_J07s_ov_08_R1.fastq.gz            Cxq_NG_J07s_ov_08_R2.fastq.gz
## Cxq_NG_J07s_ov_09_R1.fastq.gz            Cxq_NG_J07s_ov_09_R2.fastq.gz
## Cxq_NG_J14ns_ov_01_R1.fastq.gz           Cxq_NG_J14ns_ov_01_R2.fastq.gz
## Cxq_NG_J14s_ov_01_R1.fastq.gz            Cxq_NG_J14s_ov_01_R2.fastq.gz

For information, s stands for mosquitoes that have undergone salivation and ns for mosquitoes that didn’t

>ls references

wPipPel_Plasmid.fasta

grep ">" wPipPel_Plasmid.fasta

>wPipPel_GCF_000073005.1
>pWCP

AliNe is a Nextflow-type workflow, meaning we need to download the nextflow binary. No install is required since nextflow automatically downloads all dependencies using a container type system (e.g. Docker, and in our case Singularity).

curl -s https://get.nextflow.io > nextflow

We create a simple configuration file for Nextflow: nf.conf. We specify the working directory of Nextflow to a subdirectory ./tmp/ in the current directory, allowing us easy access to Nextflow’s temporary working directory and preventing any memory disk usage issues from using the cluster’s /tmp/ partitions. The dag.overwrite parameters allows to restart computing the workflow graph completely in case of any error.

# nf.conf
workDir = './tmp/'
dag.overwrite = true

Then, we created the file run_workflow.sh, for successful execution of the job on the Migale cluster with the SGE system.

We iterate the Nextflow run of AliNe over our two concatenated reference genomes: wPipPel and pWCP.

The file is as follows:

#!/bin/bash
# run_workflow.sh

# Shell to use for job execution
#$ -S /bin/bash

# Job name
#$ -N mapping

# Nom de la queue
#$ -q infinit.q

# Export all environment variables
#$ -V

# Run command from current repertory
#$ -cwd

# CPUs number
#$ -pe thread 2

# Logs
#$ -o ./workflow.out 
#$ -e ./workflow.err

# activate conda environment
conda activate anvio-8

for g in 'references/wPipPel_Plasmid.fasta'; do
        echo $g
        ./nextflow run Juke34/AliNe \
          -r v1.1.0 \
          -profile singularity \
          --reads Trimmed/ \
          --genome $g \
          --read_type short_paired \
          --aligner bowtie2 \
          --fastqc \
          --samtools_stats \
          --trimming_fastp \
          -c nf.conf
done

>ls

data    references   ./nextflow    nf.conf

We specify the read type to be short read and paired end (since we used this type of Illumina sequencing), and the aligner to be bowtie2 with default options. The version of AliNe used is 1.1.0.

Figure: Workflow process The workflow is as follows: first, we take the paired-end short reads, and send them to 1) quality control for raw data and 2) trimming after a step of score standardization. Score standardization is performed with seqkit convert, the idea of AliNe is to have a general tool for all sequencers (Illumina, Solexa, Nanopore) so this step of score standardization is used in that context. Trimming is performed with fastp. Quality control is performed with FastQC. After trimming, the trimmed reads are sent to FastQC.

The second step of the workflow is performing the alignment of the RNA reads against the reference genome. The preliminary phase of the alignment is library type guessing, which is performed with salmon quant. The goal is to know if the RNA reads are stranded or unstranded, found on the forward or reverse strand, and with reads positioned in inward, outward compared to the strand, or matching the strand. The summary of library types detected by salmon is reported here. Then alignment can be performed with any tool provided by AliNe; here in our study we chose bowtie2 as the sole aligner.

Finally, quality check is performed on the output of the alignments, with FastQC, and the results of all FastQC analyzes (before trimming, after trimming, and after alignment) are aggregated using tool MultiQC.

Once the workflow is finished, a number of folders are created, which we show here:

>ls

alignment_results   data    rna-seq-cxq   ./nextflow    nf.conf   tmp     pipeline_report   workflow.out    

The most important folder, which contains all results, is alignment_results. The pipeline visualizations are contained in pipeline_report.

>ls alignment_results

alignment  fastp  fastqc  mean_read_length  salmon_libtype  samtools_stats  seqkit_score

The .bam files generated by the bowtie2 alignment are found in alignment_results/alignment/bowtie2/:

>ls alignment_results/alignment/bowtie2/*.bam

Cxq_NG_J01_ov_10_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J01_ov_11_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J01_ov_13_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J01_ov_14_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_06_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_11_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J03_ov_12_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_06_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_08_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07ns_ov_09_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_05_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_08_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J07s_ov_09_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J14ns_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
Cxq_NG_J14s_ov_01_R1_seqkit_trim_bowtie2_sorted.bam

Here is an example of a finishing Nextflow workflow:

Figure: Example of running process of Nextflow workflow

UMI-tools deduplication

In order to remove PCR bias, our reads are furnished with Unique Molecular Identifiers (UMIs). The idea is to afford correct counting and normalization by removing duplicates from PCR amplification. This is performed with the package umi_tools : see doc.

We already have the UMI tags extracted on our data, here we perform the deduplication on our data after mapping on the .bam files.

umi_tools dedup --umi-separator=: -I [example].bam --output-stats=[deduplicated_stats].txt -S [deduplicated].bam

We use this command in a for loop on our *.bam, and obtained new bams but deduplicated and smaller in size.

In the end, we have the followign result folder:

>ls result

deduplicated_Cxq_NG_J01_ov_10_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J01_ov_11_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J01_ov_13_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J01_ov_14_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_06_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_11_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J03_ov_12_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_06_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_08_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07ns_ov_09_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_02_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_03_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_04_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_05_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_07_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_08_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J07s_ov_09_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J14ns_ov_01_R1_seqkit_trim_bowtie2_sorted.bam
deduplicated_Cxq_NG_J14s_ov_01_R1_seqkit_trim_bowtie2_sorted.bam

Final samtools treatment

In order to use the .bam files, we need to do the final samtools commands, sort and index:

for f in *.bam; do 
  samtools sort $f > ${f}_sorted.bam
done

for f in *.bam; do 
  samtools index ${f}_sorted.bam
done

Now that our .bams are ready, we can use bedtools genomecov and stringtie on the .bam to get RNA read coverage information in the next steps.