Comments to umi deduplication script #8

Merged
T. Karp merged 2 commits from s3970582/system_genetics:master into master 2021-02-23 15:24:46 +01:00
9 changed files with 361 additions and 0 deletions
Showing only changes of commit f58dfdaad6 - Show all commits

4
.gitignore vendored Normal file
View File

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

View File

@ -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}"

View File

@ -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

View File

@ -0,0 +1,54 @@
#!/bin/bash
#
# Align reads against reference genome.
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
# Store genome index in this location:.
GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
mkdir -p "${GENOME_INDEX}"
# Store the generated `Aligned.sortedByCoord.out.bam` in this dir.
ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment/"
mkdir -p "${ALIGNMENT_OUTPUT}"
# 1) Generate genome index.
#
# 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.
# - We're using gzip compressed reference data (--readFilesCommand zcat), i.e.,
# .gtf.gz and fa.gz. If not, you can remove the `zcat` flag.
# Storage location reference data (in this case on Gearshift).
REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf.gz"
FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
STAR \
--runThreadN 8 \
--runMode genomeGenerate \
--readFilesCommand zcat \
--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}"

View File

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

View File

@ -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

View File

@ -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.

View File

@ -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.