RNA-Sequencing Part 3 Generation and characterization of a novel mouse model of Becker Muscular Dystrophy with a deletion of exons 52 to 55

1. Starting with the original file (Number.Aligned.sortedByCoord.out.bam), using samtools, we sorted by name: 
 
samtools sort -n Number.Aligned.sortedByCoord.out.bam -o Number.Aligned.namesorted.bam 
 

2. We extracted the paired reads into 2 separate files for R1 and R2, and any singleton or orphaned reads into additional RS and R0 files, respectively (many of the RS and R0 files were empty and not added here due to size constraints): 
 
samtools fastq -1 Number_R1.fastq -2 Number_R2.fastq -0 Number_R0.fastq -s Number_RS.fastq 
 

3. We compressed all of the files into ‘.gz’ extension using gzip: 
 
gzip -9 Number_R1.fastq 

.bam and RS/R0 files were not added due to size constraints but were available upon request. 

Upstream workflow performed by TCAG (SickKids): 

2. RNA-Seq Library and Reference Genome Information

Type of library: stranded, paired end

Genome reference sequence: GRCm39, M31 Gencode gene models.

3. Read Pre-processing, Alignment and Obtaining Gene Counts

3.1 Read Pre-processing

The sequencing data is in FASTQ format. The quality of the data is assessed using FastQC v.0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Adaptors are trimmed using Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) v. 0.5.0. Trim Galore is running Cutadapt (https://cutadapt.readthedocs.org/en/stable/) v. 1.10. Trim Galore is run with the following parameters:

            -q 25 – the reads are trimmed from the 3' end base by base, trimming stops if the quality of the base is greater than 25;

            --clip_R1 6, --clip_R2 6 – clip the first 6 nucleotides from the 5' ends of read 1 and read 2;

            --stringency 5 – at least 5 nucleotides overlap with the Illumina primer sequence are needed for trimming;

            --length 40 – any read that is shorter than 40 nucleotides as a result of trimming is discarded;

            --paired – only pairs of reads are retained (for paired-end reads only, not for single reads).

The type of adaptor is automatically detected by screening the first 1 million sequences of the first specified file for the first 12/13 nucleotides of the standard Illumina or Nextera primers and the sequence from the start of the primer to the 3' end of the read is trimmed.

The quality of the trimmed reads is re-assessed with FastQC.

The trimmed reads are also screened for presence of rRNA and mtRNA sequences using FastQ-Screen v.0.10.0 (http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/).

To assess the read distribution, positional read duplication and to confirm the strandedness of the alignments we use the RSeQC package (http://rseqc.sourceforge.net/), v. 2.6.2. The distribution of reads across exonic, intronic and intergenic sequences is assessed by the read_distribution.py program, infer_experiment.py is used for confirming strandedness, and read_duplication.py is used to obtain the positional read duplication (percentage of reads mapping to exactly the same genomic location). Sufficient proportion of reads should map to the exonic sequences (ideally > 70-80%). Large amounts of reads mapping to intronic sequences in a poly-A mRNA library will suggest significant presence of pre-mRNA or other issues with RNA preparation. For stranded RNA-seq experiments the majority of the reads should map exclusively to one strand, same or opposite to the transcript, depending on the library preparation method. For non-stranded experiments the reads should be equally distributed to both strands.

3.2. Read Alignment

The raw trimmed reads are aligned to the reference genome using the STAR aligner, v.2.6.0c. (https://github.com/alexdobin/STAR, https://academic.oup.com/bioinformatics/article/29/1/15/272537).  The alignments are contained in the .bam files. The “.bam” together with the “.bai” files can be used for viewing of the alignments in the Integrative Genomics Viewer (IGV, http://software.broadinstitute.org/software/igv/).

3.3. Obtaining Gene Counts

The filtered STAR alignments are processed to extract raw read counts for genes using htseq-count v.0.6.1p2 (HTSeq, http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html). Assigning reads to genes by htseq-count is done in the mode “intersection_nonempty”, i.e. if a read overlaps with two overlapping genes and the overlap to gene A is greater than the overlap to gene B, the read is counted towards gene A, while if a read overlaps equally with gene A and gene B, then it is not counted towards either gene. Htseq_count does not count reads with multiple alignments to avoid introducing bias in the expression results. Only uniquely mapping reads are counted.

4. Pre-processing, Alignment and Gene Counts QC

MultiQC (https://multiqc.info/) is a reporting tool that aggregates statistics generated by bioinformatics analyses across multiple samples. MultiQC v. 1.14 was used to generate a consolidated report from FastQC screening of both untrimmed and trimmed reads, and from RSeQC,  FastQ Screen, STAR and htseq-count results. The MultiQC report is contained in MultiQC_Report_*.html file.

 5. DGE Analysis with edgeR

Differential expression was done with the edgeR R package v.3.28.1, using R v.3.6.1 (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html). The data set was filtered to retain only genes whose gene counts were >50 in at least 3 samples. This is intended to remove genes that are not expressed, or expressed at a very low level.

The method used for normalizing the data was TMM, implemented by the calcNormFactors(y) function. All samples were normalized and filtered together. The glmLRT functionality in edgeR was used for the differential expression tests, with sample group taken into account.

EdgeR Results Legend:

·      GeneID – Ensembl Gene ID;

·      Chr.Start.End - gene coordinates;

·      GeneName, GeneType, etc. – Gene attributes, derived from the genome annotation;

·      logFC - Log2 Fold Change (use this column for selection of DEGs);

·      logCPM - Log2 Counts Per Million, average for all libraries;

·      LR – Statistic calculated by the LR-Test;

·      PValue - Differential expression P value;

·      FDR – Differential expression False Discovery Rate, calculated by the Benjamini-Hochberg method (use this column for selection of DEGs);

·      (columns labeled with sample names) – Fragments Per Kilobase of transcript per Million mapped reads (FPKMs) for the given samples.