Use BAM files for htseq-count.

Reviewed-on: #8

.gitignore

@@ -0,0 +1,4 @@
+
+.DS_Store
+*.nosync
+*.RData

rnaseq/step1_fastqc/snippet.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+R1="sample1_R1.fastq.gz"
+R2="sample1_R2.fastq.gz"
+
+PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+FASTQC_OUT="${PROJECT_DIRECTORY}/step1/"
+mkdir -p "${FASTQC_OUT}"
+
+# Run FastQC on paired-end data.
+fastqc \
+    -o "${FASTQC_OUT}" \
+    "${R1}" "${R2}"

rnaseq/step2_trim/snippet.sh

@@ -0,0 +1,54 @@
+#!/bin/bash
+#
+# Reference: http://www.usadellab.org/cms/?page=trimmomatic
+
+
+module load Trimmomatic
+PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+FASTQ_OUT="${PROJECT_DIRECTORY}/step2/"
+mkdir -p "${FASTQ_OUT}"
+
+
+# Adapters can be found at
+# https://github.com/timflutre/trimmomatic/tree/master/adapters
+# But should be verified with FastQC, or in another way.
+
+
+# Trimmomatic example Paired end data.
+#
+# Flags:
+# - ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the
+#                 read.
+# - SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average
+#                  quality within the window falls below a threshold.
+# - LEADING: Cut bases off the start of a read, if below a threshold quality.
+# - TRAILING: Cut bases off the end of a read, if below a threshold quality.
+# - HEADCROP: Cut the specified number of bases from the start of the read.
+# - MINLEN: Drop the read if it is below a specified length.
+java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
+  -phred33 \
+  sample1_R1.fastq.gz \
+  sample1_R2.fastq.gz \
+  "${FASTQ_OUT}/sample1_R1_paired.fastq.gz" \
+  "${FASTQ_OUT}/sample1_R1_unpaired.fastq.gz" \
+  "${FASTQ_OUT}/sample1_R2_paired.fastq.gz" \
+  "${FASTQ_OUT}/sample1_R2_unpaired.fastq.gz" \
+  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
+  LEADING:3 \
+  TRAILING:3 \
+  SLIDINGWINDOW:4:25 \
+  HEADCROP:8 \
+  MINLEN:50
+
+
+# Example single end data.
+java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar SE \
+  -phred33 \
+  sample1.fastq.gz \
+  output.fastq.gz \
+  ILLUMINACLIP:TruSeq3-SE:2:30:10 \
+  LEADING:3 \
+  TRAILING:3 \
+  SLIDINGWINDOW:4:15 \
+  HEADCROP:8 \
+  MINLEN:50

rnaseq/step3_align/snippet.sh

@@ -0,0 +1,63 @@
+#!/bin/bash
+#
+# Align reads against reference genome.
+
+REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
+
+PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+# Store the generated `sample1_Aligned.sortedByCoord.out.bam` in this dir.
+ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment"
+mkdir -p "${ALIGNMENT_OUTPUT}"
+
+# 1) Generate genome index (optional).
+#
+# Depending on your read size, reference genome and annotation, you may need to
+# generate a new genome index. In most cases, this is not necessary and you can
+# directly use the pre-build genome index from the cluster:
+#
+# GENOME_INDEX="${REFERENCE_DATA}/index_GRCh38_gtf100_overhang100"
+#
+# and ignore the first STAR command below.
+#
+# N.B.:
+# - We're assuming a read size of 100 bp (--sjdbOverhang 100). An alternative
+#   cut-off is 150, for low-input methods. In general, refer back to the
+#   previous quality control steps if you are unsure about the size. In case of
+#   reads of varying length, the ideal value is max(ReadLength)-1.
+# - If you're using gzip compressed reference data, i.e., .gtf.gz and fa.gz,
+#   pass the `--readFilesCommand zcat` flag.
+
+# Store created genome index in this location:.
+GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
+mkdir -p "${GENOME_INDEX}"
+
+# Storage location reference data on Gearshift.
+GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
+FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa"
+
+STAR \
+    --runThreadN 8 \
+    --runMode genomeGenerate \
+    --sjdbOverhang 100 \
+    --genomeFastaFiles ${FASTA_FILE} \
+    --sjdbGTFfile ${GTF_FILE} \
+    --genomeDir ${GENOME_INDEX}
+
+
+# 2) Do the actual alignment.
+#
+# N.B.:
+# - We are assuming paired-end, gzip compressed (--readFilesCommand zcat)  FastQ
+#   files.
+
+# The compressed, paired-end, FastQ's after trimming (step 2).
+R1="${PROJECT_DIRECTORY}/step2/sample1_R1_paired.fastq.gz"
+R2="${PROJECT_DIRECTORY}/step2/sample1_R2_paired.fastq.gz"
+
+STAR \
+    --runThreadN 8 \
+    --readFilesCommand zcat \
+    --readFilesIn "${R1}" "${R2}" \
+    --outSAMtype BAM SortedByCoordinate \
+    --genomeDir ${GENOME_INDEX} \
+    --outFileNamePrefix "${ALIGNMENT_OUTPUT}/sample1_"

rnaseq/step4_multiqc/snippet.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+module load multiqc
+
+PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+multiqc "${PROJECT_DIRECTORY}"

rnaseq/step5_count/snippet.sh

@@ -0,0 +1,28 @@
+#!/bin/bash
+PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
+COUNT_OUTPUT="${PROJECT_DIRECTORY}/step5"
+mkdir -p "${COUNT_OUTPUT}"
+
+# Storage location of annotation on Gearshift.
+REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
+GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
+
+# Where our alignment file was stored.
+BAM="${PROJECT_DIRECTORY}/step3/alignment/sample1_Aligned.sortedByCoord.out.bam"
+
+# Compute counts using htseq-count.
+#
+# N.B.:
+# - If you are processing multiple files, consider using the `--nprocesses` flag
+#   to distribute computation of the files to different cores.
+# - The BAM file must be position sorted. If you used STAR with the
+#   `SortedByCoordinate` option you should be okay. If not, sort your BAM file
+#   using `samtools sort`.
+# - By default, strand aware library preparation is assumed. If not, specify the
+#   `--stranded` flag.
+htseq-count \
+    --order pos \
+    -f bam \
+    ${BAM} \
+    ${GTF_FILE} \
+    > ${COUNT_OUTPUT}/counts.txt

rnaseq/step6_overall_QC/01 - Principle Component Analysis.R

@@ -0,0 +1,81 @@
+# Principle Component Analysis
+# Normalized with limma::voom
+library(limma)
+library(tidyverse)
+
+
+# expression.data. Has columns:
+# - Gene (gene identifier)
+# - [Sample identifiers]
+expression.data <- "mrna_counts_table.csv"
+# results.dir. The directory of where to put the resulting tables (and later your plots.)
+results.dir <- "results"
+
+
+# PCA variables
+do.center = TRUE
+do.scale = FALSE
+
+
+# The analysis
+# We use prcomp to calculate the PCAs. Afterwards you should plot the results.
+norm.expr.data <- expression.data %>%
+	tibble::column_to_rownames("Gene")
+norm.expr.data <- norm.expr.data[rowSums(norm.expr.data) >= 10,] %>%
+	limma::voom() %>%
+	as.matrix()
+
+# Principle Component analysis
+results.dir.pca <- file.path(results.dir, "principle.components")
+dir.create(results.dir.pca, recursive=TRUE)
+
+norm.expr.data.pcs <- norm.expr.data %>%
+  t() %>%
+  stats::prcomp(
+  	center = do.center,
+  	scale. = do.scale
+  )
+
+# Write summary of PCAs to files
+pcs.summary <- summary(norm.expr.data.pcs)
+pcs.summary$importance %>%
+	t() %>%
+	as.data.frame() %>%
+	tibble::rownames_to_column("PC.name") %>%
+	readr::write_csv(
+		file.path(results.dir.pca, "importance.csv")
+	)
+
+pcs.summary$x %>%
+	t() %>%
+	as.data.frame() %>%
+	tibble::rownames_to_column("ensembl.id") %>%
+	readr::write_csv(
+		file.path(results.dir.pca, "values.csv")
+	)
+
+pcs.summary$rotation %>%
+	t() %>%
+	as.data.frame() %>%
+	tibble::rownames_to_column("sample.id") %>%
+	readr::write_csv(
+		file.path(results.dir.pca, "rotation.csv")
+	)
+
+data.frame(
+		rownames = names(pcs.summary$center),
+		center = pcs.summary$center,
+		scale = pcs.summary$scale
+	) %>%
+	readr::write_csv(
+		file.path(results.dir.pca, "rest.csv")
+	)
+
+# Not saved: pcs.summary$sdev,
+
+
+# Next thing to do:
+# - (Optional) scree plot - to determine the optimal cutoff for PCA inclusion based on explaination of variance
+# - (Optional) eigencorplot - to correlate PCAs to clinical variables so that you know which PCA to include for which analysis
+# - (Optional) pairsplot - plot multiple PCAs against each other in a single figure
+# - Plot the first couple of PCAs against each other

rnaseq/step6_overall_QC/02 - Gender Check.R

@@ -0,0 +1,112 @@
+# Gender QC
+# Normalized with limma::voom
+library(GSVA)
+library(limma)
+library(edgeR)
+library(tidyverse)
+
+
+# expression.data. Has columns:
+# - Gene (gene identifier)
+# - [Sample identifiers]
+expression.data <- "mrna_counts_table.csv"
+# master.Table. Has columns:
+# - GenomeScan_ID
+# - gender, levels = c("male", "female")
+# - age
+# - factor(smoking.status, levels = c("Ex-smoker", "Current smoker"))
+master.Table <- "patient_table.csv"
+# results.dir. The directory of where to put the resulting tables (and later your plots.)
+results.dir <- "results"
+
+
+# The analysis
+# We first do a differential expression analysis on gender using EdgeR.
+# Afterwards you should plot these results.
+x.genes <- gene.data %>%
+    dplyr::filter(chromosome_name == "X") %>%
+    dplyr::pull(ensembl_gene_id)
+
+y.genes <- gene.data %>%
+    dplyr::filter(chromosome_name == "Y") %>%
+    dplyr::pull(ensembl_gene_id)
+
+# Gender QC
+results.dir.gender <- file.path(results.dir, "gender.check")
+dir.create(results.dir.gender, recursive=TRUE)
+
+# Differential Expression on Gender
+gender.qc.patients <- master.Table %>%
+    dplyr::filter(
+        !is.na(GenomeScan_ID)
+    ) %>%
+    dplyr::mutate(
+        gender = factor(gender, levels = c("male", "female")),
+        age = as.numeric(age),
+        smoking.status = factor(smoking.status, levels = c("Ex-smoker", "Current smoker"))
+    ) %>%
+    dplyr::filter(
+        !is.na(gender)
+    )
+gender.qc.sample.order <- gender.qc.patients %>%
+    dplyr::pull(GenomeScan_ID)
+
+gender.qc.expression.data <- expression.data %>%
+    tibble::column_to_rownames("Gene") %>%
+    select.columns.in.order(gender.qc.sample.order) %>%
+    as.matrix()
+
+design <- model.matrix( ~0 + gender, data = gender.qc.patients)
+
+DGEL <- edgeR::DGEList(gender.qc.expression.data)
+keep <- edgeR::filterByExpr(DGEL)
+keep[names(keep) %in% x.genes] <- TRUE
+keep[names(keep) %in% y.genes] <- TRUE
+DGEL <- DGEL[keep, , keep.lib.sizes=FALSE]
+DGEL <- edgeR::calcNormFactors(DGEL, method = "TMM")
+
+DGEL <- edgeR::estimateDisp(DGEL, design)
+fit <- edgeR::glmQLFit(DGEL,design)
+
+contrasts <- limma::makeContrasts(
+    gender = gendermale - genderfemale,
+    levels = design
+)
+qlf   <- edgeR::glmQLFTest(fit, contrast = contrasts[,"gender"])
+
+gender.qc.results <- edgeR::topTags(
+        qlf,
+        n=nrow(DGEL)
+    )$table %>%
+    tibble::rownames_to_column("ensembl.id") %>%
+    dplyr::left_join(
+        y = gene.data,
+        by = c("ensembl.id" = "ensembl_gene_id")
+    ) %>%
+    readr::write_csv(
+        file.path(results.dir.gender, "differential.expression.on.gender.csv")
+    )
+
+# Plotting of gender expression
+results.dir.gender.plot <- file.path(results.dir.gender, "img")
+dir.create(results.dir.gender.plot, recursive=TRUE)
+
+gender.qc.genes.to.plot <- gender.qc.results %>%
+    dplyr::arrange(PValue) %>%
+    dplyr::group_by(chromosome_name) %>%
+    dplyr::filter(
+        (
+            chromosome_name %in% c("X", "Y") &
+            FDR < 0.05 &
+            dplyr::row_number() <= 5
+        )
+    )
+
+
+# Next thing to do:
+# - Plot the normalized expressino values for the genes in gender.qc.genes.to.plot in a boxplot, split and colored by gender.
+# - (Optional) Do a GSVA with as genesets the genes found in gender.qc.genes.to.plot. Plot the boxplots as per the previous point.
+# - (Optional) Plot the number of Y-chromosome reads devided by the number of X chromosome reads in a boxplot as per the first point.
+#     x=$(samtools view -q 20 -f 2 $bam_file X | wc -l)
+#     y=$(samtools view -q 20 -f 2 $bam_file Y | wc -l)
+# - (Optional) Plot the number of Y-chromosome SNPs devided by the number of X chromosome SNPs in a boxplot as per the first point.

rnaseq/step6_overall_QC/03 - Sample Counts.R

@@ -0,0 +1,39 @@
+# Total counts per sample
+# Normalized with limma::voom
+library(limma)
+library(tidyverse)
+
+
+# expression.data. Has columns:
+# - Gene (gene identifier)
+# - [Sample identifiers]
+expression.data <- "mrna_counts_table.csv"
+# master.Table. Has columns:
+# - GenomeScan_ID
+# - gender, levels = c("male", "female")
+# - age
+# - factor(smoking.status, levels = c("Ex-smoker", "Current smoker"))
+master.Table <- "patient_table.csv"
+# results.dir. The directory of where to put the resulting tables (and later your plots.)
+results.dir <- "results"
+
+
+# The analysis
+# We calculate the number of mapped reads per sample.
+total.count.per.sample <- expression.data %>%
+	tibble::column_to_rownames("Gene") %>%
+	colSums()
+
+data.frame(
+		sample = names(total.count.per.sample),
+		counts = as.numeric(total.count.per.sample)
+	) %>%
+	readr::write_csv(file.path(results.dir, "total.counts.per.sample.csv"))
+
+
+# Next thing to do:
+# - Check the number of reads per sample in total.counts.per.sample.csv
+# - Plot the reads distribution (all reads) per sample in a boxplot.
+# - (Optional) Calculate the number of unmapped, multimapped, unique mapped to
+#   feature and unique mapped to no feature and plot these in a stacked bar graph.
+

rnaseq/umi/07_remove_dups.pl

@@ -0,0 +1,61 @@
+#!/usr/bin/perl -w
+use strict;
+use Parallel::ForkManager;
+# this script creats one file with UMI unique reads and one with UMI duplicated reads
+# as input you need aligned sorted by coordinate bam file
+my @torun = ();
+foreach my $file ( <*Aligned.sortedByCoord.out.bam> ) {
+    push @torun, $file;
+}
+
+my $pm = Parallel::ForkManager->new( 12 );
+foreach my $file ( @torun ) {
+    my $sample = $file;
+    $sample =~ s/Aligned\.sortedByCoord\.out\.bam$//;
+    next if -s $sample.'_uniq.bam';
+    warn "Parsing $sample\n";
+    $pm->start and next;
+    my %seen = ();
+    my %duplicates = ();
+    my ( $uniqs, $dups ) = ( 0, 0 );
+    open F, 'samtools view -h '.$file.' |';
+    open F1, '| samtools view -bS - >'.$sample.'_uniq.bam';
+    open F2, '| samtools view -bS - >'.$sample.'_dups.bam';
+    while ( <F> ) {
+        if ( m/^\@/ ) {
+            print F1;
+            print F2;
+            next;
+        }
+        my ( $id, $flag, $chr, $pos, $mapq, $cigar, $chr2, $pos2, $tlen ) = split /\t/;
+        next if $flag & 256 or $flag & 512 or $flag & 1024; #skip if the read is not primary alignment/read fails platform/vendor quality checks/read is PCR or optical duplicate
+#        foreach ( 256, 512, 1024 ) { $flag-=$_ if $flag&$_ }
+        my ( $bc ) = $id =~ m/\:([GATCN\d]+)$/; #extract UMI barcode
+        my $uniq = join( ':', $chr, $pos, $flag, $tlen, $bc );
+        my $pos_ = $pos-1;
+        while ( $cigar =~ m/(\d+)([SHMDIN=])/g ) {
+            $pos_+=$1 if $2 eq 'M' or $2 eq '=' or $2 eq 'D' or $2 eq 'N'; # find position for minus strand
+        }
+        my $uniq2 = join( ':', $chr, $pos_, $flag, $tlen, $bc );
+        if ( exists($duplicates{$id}) or # already marked as duplicate
+            ( not($flag&16) and ++$seen{ $uniq } > 1 ) or # plus strand
+            ( $flag&16 and ++$seen{ $uniq2 } > 1 ) #minus strand
+            ) {
+            print F2;
+            $dups++;
+            $duplicates{$id} = 1;
+        }
+        else {
+            print F1;
+            $uniqs++;
+        }
+    }
+    close F;
+    close F1;
+    close F2;
+    warn "Unique: $uniqs (", 100*$uniqs/($uniqs+$dups), ")\n";
+    system( 'samtools', 'index', $sample.'_uniq.bam' );
+#    last;
+    $pm->finish;
+}
+$pm->wait_all_children;