Tutorial of TSS differential analysis
Differential analysis using refTSS
This use-case provides step by step workflow to achieve gold standard expression analysis using refTSS, and obtains TSS expression table, input TSSID list for refTSS-LD, quality assessment of TSS sequencing and statistical testing results. These scripts are based on the scripts written in Morioka MS et al..
For hands on tutorial, refTSS provides an example data-set published by Yasukawa et al. at use-case file folder:
Software Requirements
- NCBI tools (for fastq file download)
- bigWigAverageOverBed, bedGraphToBigWig in Jim Kent Utility tools
- fastqc program
- cutadapt
- STAR
- R (version 4.2.2 or later)
- bedtools (Confirmed to work with version 2.28 and version 2.29)
1. FASTQ preparation
1.1 Download genome reference files
GRCh38 human genome assembly file is downloaded from GENCODE. After downloading genome fasta file, `.genome` file which required in several utility tools is generated by `faidx` function in NCBI tools.
# .genome file generation pip install --upgrade pip pip install pyfaidx gzip -d GRCh38.p14.genome.fa.gz faidx GRCh38.p14.genome.fa -i chromsizes >GRCh38.p14.genome
1.2 FASTQ file preparation
Downloading an example data set (CAGE data) from SRA.
fasterq-dump DRR159241 -o mutA_rep1.fq -e 12 -p fasterq-dump DRR159242 -o mutA_rep2.fq -e 12 -p fasterq-dump DRR159243 -o mutA_rep3.fq -e 12 -p fasterq-dump DRR159244 -o wt_rep1.fq -e 12 -p fasterq-dump DRR159245 -o wt_rep2.fq -e 12 -p fasterq-dump DRR159246 -o wt_rep3.fq -e 12 -p
Sequencing quality check by fastqc program
fastqc -o qc_mutA_rep1 mutA_rep1.fq fastqc -o qc_mutA_rep2 mutA_rep2.fq fastqc -o qc_mutA_rep3 mutA_rep3.fq fastqc -o qc_wt_rep1 wt_rep1.fq fastqc -o qc_wt_rep2 wt_rep2.fq fastqc -o qc_wt_rep3 wt_rep3.fq
1.3 [Optional] Adaptor trimming step (e.g., two color illumina chemistries)
One of the major sequencing techs are both NextSeq and the NovaSeq which uses two-colored chemisty. In this case, 3'-end of `G` was derived from end of sequencing with good sequence quality. cutadapt with `--nextseq-trim=20` option will work as regular quality trimming (where one would use `-q 20` instead), except that the qualities of 3'-end `G` bases are ignored.
# -m: minimum output read length # -a: adaptor trimming # Trueseq read1 adaptor sequence:AGATCGGAAGAGCACACGTCTGAACTCCAGTCA # Nextra adaptor sequence: CTGTCTCTTATACACATCT # --nextseq-trim=20 cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o mutA_rep1.tm.fq mutA_rep1.fq cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o mutA_rep2.tm.fq mutA_rep2.fq cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o mutA_rep3.tm.fq mutA_rep3.fq cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o wt_rep1.tm.fq wt_rep1.fq cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o wt_rep2.tm.fq wt_rep2.fq cutadapt -m 30 --nextseq-trim=20 -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA -o wt_rep3.tm.fq wt_rep3.fq
2. Mapping TSS-seq reads on reference genome
This example uses STAR program for mapping read on reference genome.
# STAR genome index generation
STAR --runMode genomeGenerate \
--genomeDir GRCh38p14_star_index \
--runThreadN 16 \
--genomeFastaFiles GRCh38.p14.genome.fa \
--sjdbGTFfile gencode.v46.annotation.gtf
# STAR Mapping with example options
GTF=/home/masaki.morioka/genome/gencode.v46.annotation.gtf
SGENOME=/home/masaki.m/genome/GRCh38p14_star_index
_naming=mutA_rep1
STAR --runThreadN 12 --genomeDir $SGENOME \
--readFilesIn ${_naming}.tm.fq \
--outSAMtype BAM SortedByCoordinate \
--outTmpKeep None \
--outFileNamePrefix ${_naming}
3. CTSS file generation and quality assessment of TSS sequencing
3.1 Download and running the script for TSS tag count table
Download the our custom script for counting TSS tags and the reference of canonical TSS coordinates (GENCODE).
# Download scripts from Github # git clone https://github.com/suimye/cage_tutorial.git chmod 755 ./cage_tutorial/* # Download canonical promoter coordinates (e.g., GENCODE) wget https:://refTSS.riken.jp/tutorial/gencode.v40.promoter300-100.bed.gz gzip -d gencode.v40.promoter300-100.bed.gz # Making a directory for tag count result mkdir gencode_count cd gencode_count/
Running the script for tag counting
# Example: bash cage.counting.pipeline.b0.01.sh {sample.bam} {mapping quality (integer)} {referece.bed} {genome size file} =
sh ../cage_tutorial/cage.counting.pipeline.b0.01.sh ../mutA_rep1Aligned.sortedByCoord.out.bam 20 gencode.v40.promoter300-100.bed ../GRCh38.p14.genome
-----
Following six tag count files are generated.
- wt_rep1Aligned.sortedByCoord.out.240925.ctss.txt - wt_rep2Aligned.sortedByCoord.out.240925.ctss.txt - wt_rep3Aligned.sortedByCoord.out.240925.ctss.txt - mutA_rep1Aligned.sortedByCoord.out.240925.ctss.txt - mutA_rep2Aligned.sortedByCoord.out.240925.ctss.txt - mutA_rep3Aligned.sortedByCoord.out.240925.ctss.txt
3.2 TSS Quality Evaluation and Create Expression Tables Using R
TSS sequencing technologies are expected to enrich TSS information from sequencing reads. Yield of tags should be overlapped with canonical TSS annotation such as 5 prime end of GENCODE mRNA annotation.
The script `promoter_mapping_rate.R` reads the files ending with `ctss.txt` generated in the previous process and creates an expression data table. Be careful not to have any unnecessary `ctss.txt` files in your working folder as it may cause errors. After the args option, enter the labels for the comparison groups (in this case, WT and Mutant) separated by spaces. These input labels will use for regular expression of file names. If you specify "none" for the `--args` option, it is used when you want to simply create an expression data table without considering replicate experiments.
= Generating normalized data for two groups using RLE normalization in edgeR = R --slave --vanilla --args wt mutA < ../cage_tutorial/promoter_mapping_rate.b0.2.R = simply generating the TSS tag count matrix = R --slave --vanilla --args none < ../cage_tutorial/promoter_mapping_rate.b0.2.R
Here is the example output `TSS.qcbarplot.pdf` from result of test data set.
4. Differential expression analysis of TSS sequencing using refTSS
TSS tags will be useful for gene expression analysis. Calculating TSS tag counts on refTSS will provide TSS-wide comprehensive view of expression changes in the genome. This process use shell script and R.
Making directory for this analysis
# make new directory cd ../ mkdir reftss_count cd reftss_count # download refTSS human corrdinate file wget https://reftss.riken.jp/datafiles/4.1/human/refTSS_v4.1_human_coordinate.hg38.bed.txt.gz gzip -d refTSS_v4.1_human_coordinate.hg38.bed.txt.gz
4.1 Calculation of TSS tag counts on refTSS
# Example: bash cage.counting.pipeline.b0.02.sh {sample.bam} {mapping quality (integer)} {referece.bed} {genome size file}
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../mutA_rep1Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../mutA_rep2Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../mutA_rep3Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../wt_rep1Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../wt_rep2Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
sh ../cage_tutorial/cage.counting.pipeline.b0.02.sh ../wt_rep3Aligned.sortedByCoord.out.bam 20 refTSS_v4.1_human_coordinate.hg38.bed.txt ../GRCh38.p14.genome
4.2 Making TSS expression matrix in R environment
R --slave --vanilla --args wt mutA < ../cage_tutorial/promoter_mapping_rate.b0.2.R
After execution of previous command, R data file was generated as `cage_count{date}.RData`
4.3 Statistical testing
Previous process generate a R data file in working directory. Loading R data file and calculating statistics on R console.
library(edgeR)
load("cage_count{date}.RData")
# R object list
> ls()
[1] "args" "cname" "d" "df" "dgeL"
[6] "fl" "flag" "gp" "gp1" "gp2"
[11] "i" "log2cpm" "read_d" "rname" "size"
[16] "time_stamp"
# Statistical testing by exact test
et_res <- exactTest(dgeL)
# Extraction of top 100 TSSs by statistics (p-value, FDR)
detss <- topTags(et_res,n=100,adjust.method="BH",sort.by="Pvalue",p.value=0.0001)
# Dividing up- and down-regulated TSS list
de_up <- subset(detss$table,logFC>0)
de_down <- subset(detss$table,logFC<0)
4.4 Visualization of differential expression by MA-plot
One of the standard approaches on visualization of expression value is MA-plot showing relationship between average expression and fold difference in each genes (or TSSs). Here is an example of MA-plot using test data.
pdf("maplot.pdf")
plot(et_res$table$logCPM,et_res$table$logFC,cex=0.3,pch=16,col="darkgrey",xlab="Average log-expression",ylab="Expression log-ratio (MutA vs WT)")
points(de_up$logCPM,de_up$logFC,cex=0.7,pch=16,col="deeppink")
points(de_down$logCPM,de_down$logFC,cex=0.7,pch=16,col="steelblue")
abline(h=0,col="tan4")
dev.off()
4.5 GWAS-LD enrichment analysis using differentially expressed TSS information
One of the standard approaches in downstream analysis of gene expression is enrichment analysis that can be useful for evaluating biological function under condition. To prepare input for GWAS-LD enrichment analysis on refTSS, TSS ID lists were restored in working directory. Algorithm of GWAS-LD enrichment analysis is described in here.
# This process removes double quotation and line numbers from rowname vectors.
# upregulated list
sink("up_det.list.txt")
cat(rownames(de_up), sep = "\n")
sink()
# downregulated list
sink("down_det.list.txt")
cat(rownames(de_down), sep = "\n")
sink()
- Access to https://reftss.riken.jp/reftss/GWAS-LD_enrichment - input text from up_det.list.txt and down_det.lsit.txt - Push run botton