2020-10-24 RNAseq 从fq开始分析全流程
2020-10-24 本文已影响0人
程凉皮儿
分析流程
image.png
1.上传四个样本原始就文件到服务器
参考https://www.jianshu.com/p/a84cd44bac67从原始数据开始分析
参考https://github.com/twbattaglia/RNAseq-workflow
Step 1. Analysing Sequence Quality with FastQC
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Description
"FastQC aims to provide a simple way to do some quality control checks on raw sequence data coming from high throughput sequencing pipelines. It provides a modular set of analyses which you can use to give a quick impression of whether your data has any problems of which you should be aware before doing any further analysis."
The first step before processing any samples is to analyze the quality of the data. Within the fastq file is quality information that refers to the accuracy (% confidence) of each base call. FastQC looks at different aspects of the sample sequences to determine any irregularies or features that make affect your results (adapter contamination, sequence duplication levels, etc.)
Installation
conda install -c bioconda fastqc --yes
Command
# Help
fastqc -h
# Run FastQC
fastqc \
-o results/1_initial_qc/ \
--noextract \
input/sample.fastq
Output
── results/1_initial_qc/
└── sample_fastqc.html <- HTML file of FastQC fquality analysis figures
└── sample_fastqc.zip <- FastQC report data
2.质控
fastqc生成质控报告，multiqc将各个样本的质控报告整合为一个。
(rna) cheng@super-Super-Server:~/project$ ls
NALM6-ctrl_combined_R1.filtered.1.fq.gz
NALM6-ctrl_combined_R1.filtered.2.fq.gz
NALM6-treat_combined_R1.filtered.1.fq.gz
NALM6-treat_combined_R1.filtered.2.fq.gz
RS411-ctrl_combined_R1.filtered.1.fq.gz
RS411-ctrl_combined_R1.filtered.2.fq.gz
RS411-treat_combined_R1.filtered.1.fq.gz
RS411-treat_combined_R1.filtered.2.fq.gz
(rna) cheng@super-Super-Server:~/project$ ls *gz | xargs fastqc -t 2
第一步fastqc生成质控报告，-t设为2时同时分析2个原始fq数据，等待约30min完成，因为本服务器没安装multiqc，就暂时不做整合。
Step2: filter the bad quality reads and remove adaptors.
- 运行如下代码，得到名为config的文件，包含两列数据
mkdir $wkd/clean
cd $wkd/clean
ls /home/data/cheng/project/raw/*1.fq.gz >1
ls /home/data/cheng/project/raw/*2.fq.gz >2
paste 1 2 > config
(rna) cheng@super-Super-Server:~/project$ pwd
/home/cheng/project
(rna) cheng@super-Super-Server:~/project$ mkdir -p clean
(rna) cheng@super-Super-Server:~/project$ cd clean/
(rna) cheng@super-Super-Server:~/project/clean$ ls /home/data/cheng/project/raw/*.1.fq.gz >1
(rna) cheng@super-Super-Server:~/project/clean$ ls /home/data/cheng/project/raw/*.2.fq.gz >2
(rna) cheng@super-Super-Server:~/project/clean$ paste 1 2 > config
(rna) cheng@super-Super-Server:~/project/clean$ cat config
/home/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.1.fq.gz /home/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.2.fq.gz
/home/cheng/project/raw/NALM6-treat_combined_R1.filtered.1.fq.gz /home/cheng/project/raw/NALM6-treat_combined_R1.filtered.2.fq.gz
/home/cheng/project/raw/RS411-ctrl_combined_R1.filtered.1.fq.gz /home/cheng/project/raw/RS411-ctrl_combined_R1.filtered.2.fq.gz
/home/cheng/project/raw/RS411-treat_combined_R1.filtered.1.fq.gz /home/cheng/project/raw/RS411-treat_combined_R1.filtered.2.fq.gz
(rna) cheng@super-Super-Server:~/project/clean$ vi qc.sh
(rna) cheng@super-Super-Server:~/project/clean$ bash qc.sh config
- 打开文件 qc.sh ，并且写入如下内容
trim_galore，用于去除低质量和接头数据
source activate rna
bin_trim_galore=trim_galore
dir='/home/data/cheng/project/clean'
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
$bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
source deactivate
- 运行qc.sh
bash qc.sh config #config是传递进去的参数
- 结果显示如下
RUN STATISTICS FOR INPUT FILE: /home/cheng/project/raw/RS411-ctrl_combined_R1.filtered.1.fq.gz
=============================================
31046571 sequences processed in total
The length threshold of paired-end sequences gets evaluated later on (in the validation step)
Writing report to '/home/cheng/project/clean/RS411-ctrl_combined_R1.filtered.2.fq.gz_trimming_report.txt'
SUMMARISING RUN PARAMETERS
==========================
Input filename: /home/cheng/project/raw/RS411-ctrl_combined_R1.filtered.2.fq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.4_dev
Cutadapt version: 2.6
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
Output file(s) will be GZIP compressed
Cutadapt seems to be fairly up-to-date (version 2.6). Setting -j -j 1
Writing final adapter and quality trimmed output to RS411-ctrl_combined_R1.filtered.2_trimmed.fq.gz
>>> Now performing quality (cutoff '-q 25') and adapter trimming in a single pass for the adapter sequence: 'AGATCGGAAGAGC' from file /home/cheng/project/raw/RS411-ctrl_combined_R1.filtered.2.fq.gz <<<
Connection to 192.168.1.6 closed by remote host.
Connection to 192.168.1.6 closed.
我的程序被中断了，因为内存不足可能
管理员给我发一个提醒：
(rna) cheng@super-Super-Server:~$ cat A_Message_from_Admin
your programme almost used all the disk, so I killed it.
Please move all your data to /home/data, this directory is on the other disk and it is big enough to run your programme.
After you read this message, you can delete it by yourself. I changed the owner to you.
更换工作目录：
(rna) cheng@super-Super-Server:~$ mv project/ /home/data/
(rna) cheng@super-Super-Server:~$ cd /home/data/
(rna) cheng@super-Super-Server:/home/data$ mv project/ cheng/
(rna) cheng@super-Super-Server:/home/data$ cd cheng
(rna) cheng@super-Super-Server:/home/data/cheng$ ls
project
查看目录下软件信息
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/
genome gtf hisat index trnadb
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/biosoft/
cellranger-4.0.0 JAGS-4.3.0 ncbi-blast-2.10.1+ refdata-gex-GRCh38-2020-A.tar.gz tophat-2.1.1.Linux_x86_64.tar.gz
cellranger-4.0.0.tar.tar JAGS-4.3.0.tar.gz ncbi-blast-2.10.1+-x64-linux.tar.gz RNA-SeQC wget.out
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/biosoft/RNA-SeQC/
gencode.v7.annotation_goodContig.gtf.gz gencode.v7.gc.txt Homo_sapiens_assembly19.fasta.gz rRNA.tar.gz ThousandReads.bam
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/genome/
GRCh38_reference_genome hg19 hg38 M25 mm10
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/gtf/
ensembl gencode
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/hisat/
grcm38.tar.gz hg19.tar.gz hg38.tar.gz mm10.tar.gz
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/index/
bowtie bwa hisat tophat
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/trnadb/
hg38 mm10
本服务器可以用比对软件较多：
(rna) cheng@super-Super-Server:/home/data/cheng$ ls /home/data/reference/index/
bowtie bwa hisat tophat
先完成去去接头:
#确定工作目录
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ pwd
/home/data/cheng/project/clean
#构建 config
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ ls /home/data/cheng/project/raw/*1.fq.gz >1
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ ls /home/data/cheng/project/raw/*2.fq.gz >2
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ paste 1 2 > config
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ ls
1 2 config
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ cat config
/home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.2.fq.gz
/home/data/cheng/project/raw/NALM6-treat_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/NALM6-treat_combined_R1.filtered.2.fq.gz
/home/data/cheng/project/raw/RS411-ctrl_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/RS411-ctrl_combined_R1.filtered.2.fq.gz
/home/data/cheng/project/raw/RS411-treat_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/RS411-treat_combined_R1.filtered.2.fq.gz
#编写质控脚本
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ vi qc.sh
#将任务挂到后台运行
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ nohup bash qc.sh config &
[1] 16284
step4: alignment
star, hisat2, bowtie2, tophat, bwa, subread都是可以用于比到的软件
可参考https://www.jianshu.com/p/a84cd44bac67
- 先运行一个样本，测试一下
mkdir $wkd/test
cd $wkd/test
source activate rna
ls $wkd/clean/*gz |while read id;do (zcat ${id}|head -1000> $(basename ${id} ".gz"));done
id=SRR1039508
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq -2 ${id}_2_val_2.fq -S ${id}.hisat.sam
subjunc -T 5 -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq -R ${id}_2_val_2.fq -o ${id}.subjunc.sam
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38 -1 ${id}_1_val_1.fq -2 ${id}_2_val_2.fq -S ${id}.bowtie.sam
bwa mem -t 5 -M /public/reference/index/bwa/hg38 ${id}_1_val_1.fq ${id}_2_val_2.fq > ${id}.bwa.sam
事实上，对RNA-seq数据来说，不要使用bwa和bowtie这样的软件，它需要的是能进行跨越内含子比对的工具。所以后面就只先择hisat2，STAR进行下一步，但是本服务器中hisat2的索引文件不完整，最后只用STAR进行了比对
- 批量比对代码
cd $wkd/clean
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_combined_R1.filtered.1_val_1.fq.gz ${id}_combined_R1.filtered.2_val_2.fq.gz
hisat2 -p 10 -x /home/data/reference/index/hisat/hg38/genome -1 ${id}_combined_R1.filtered.1_val_1.fq.gz -2 ${id}_combined_R1.filtered.2_val_2.fq.gz -S ${id}.hisat.sam
done
Note the two inputs for this command are the genome located in the (genome/ folder) and the annotation file located in the (annotation/ folder)
# This can take up to 30 minutes to complete
STAR \
--runMode genomeGenerate \
--genomeDir star_index \
--genomeFastaFiles /home/data/reference/genome/hg38/* \
--sjdbGTFfile /home/data/reference/gtf/ensembl/* \
--runThreadN 4
# Help
STAR -h
# Run STAR (~10min)
STAR \
--genomeDir star_index \
--readFilesIn filtered/sample_filtered.fq \
--runThreadN 4 \
--outSAMtype BAM SortedByCoordinate \
--quantMode GeneCounts
# Move the BAM file into the correct folder
mv -v results/4_aligned_sequences/sampleAligned.sortedByCoord.out.bam results/4_aligned_sequences/aligned_bam/
# Move the logs into the correct folder
mv -v results/4_aligned_sequences/${BN}Log.final.out results/4_aligned_sequences/aligned_logs/
mv -v results/4_aligned_sequences/sample*Log.out results/4_aligned_sequences/aligned_logs/
先构建STAR索引
Generating Indexes
Similar to the SortMeRNA step, we must first generate an index of the genome we want to align to, so that there tools can efficently map over millions of sequences. The star_index folder will be the location that we will keep the files necessary to run STAR and due to the nature of the program, it can take up to 30GB of space. This step only needs to be run once and can be used for any subsequent RNAseq alignment analyses.
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ STAR
Usage: STAR [options]... --genomeDir /path/to/genome/index/ --readFilesIn R1.fq R2.fq
Spliced Transcripts Alignment to a Reference (c) Alexander Dobin, 2009-2020
STAR version=2.7.5a
STAR compilation time,server,dir=Tue Jun 16 12:17:16 EDT 2020 vega:/home/dobin/data/STAR/STARcode/STAR.master/source
For more details see:
<https://github.com/alexdobin/STAR>
<https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf>
To list all parameters, run STAR --help
(rna) cheng@super-Super-Server:/home/data/cheng/project/clean$ cd /home/data/cheng
(rna) cheng@super-Super-Server:/home/data/cheng$ ls
project
(rna) cheng@super-Super-Server:/home/data/cheng$ nohup STAR \
> --runMode genomeGenerate \
> --genomeDir star_index \
> --genomeFastaFiles /home/data/reference/genome/* \
> --sjdbGTFfile /home/data/reference/gtf/* \
> --runThreadN 4 &
[2] 18522
(rna) cheng@super-Super-Server:/home/data/cheng$ nohup: ignoring input and appending output to 'nohup.out'
[2]+ Exit 104 nohup STAR --runMode genomeGenerate --genomeDir star_index --genomeFastaFiles /home/data/reference/genome/* --sjdbGTFfile /home/data/reference/gtf/* --runThreadN 4
(rna) cheng@super-Super-Server:/home/data/cheng$ ls
nohup.out project star_index _STARtmp
参考https://blog.csdn.net/yssxswl/article/details/105703869完成索引构建
STAR --runThreadN 20 --runMode genomeGenerate \
--genomeDir /home/data/cheng/star_index \
--genomeFastaFiles /home/data/reference/genome/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa \
--sjdbGTFfile /home/data/reference/gtf/gencode.v35.annotation.gtf \
–runThreadN 是指构建是使用的线程数，在没有其他数据在跑的情况下，可以满线程跑
–runMode genomeGenerate 让STAR执行基因组索引的生成工作
–genomeDir 构建好的参考基因组存放的位置，最好是单独建立的一个文件夹
–genomeFastaFiles 参考基因组序列文件
–sjdbGTFfile 基因注释文件
比对
我们来看一下比对的代码
STAR
--runThreadN 20 \
--genomeDir /home/data/cheng/star_index \
--readFilesCommand gunzip -c \
--readFilesIn /home/data/cheng/ /home/data/cheng/ \
--outSAMtype BAM SortedByCoordinate \
--outFileNamePrefix N052611_Alb \
–runThreadN 运行的线程数，根据自己的服务器合理选择
–genomeDir 构建的参考基因组位置
–readFilesCommand 对于gz压缩的文件，我们可以在后面添加 gunzip -c
–readFilesIn 输入文件的位置，对于双末端测序文件，用空格分隔开就行了
–outSAMtype 默认输出的是sam文件，我们这里的BAM SortedByCoordinate是让他输出为ban文件，并排序
–outFileNamePrefix 表示的是输出文件的位置和前缀
然后就是输出文件的问题，输出的文件不止一个，包含了比对过程中的一些信息
- Aligned.out.sam或者Aligned.out.bam
它指的就是我们的比对结果 - Log.progress.out
它是每分钟记录一次的对比情况 - Log.out
它记录了STAR程序在运行中的各种情况，当我们的结果出现异常时，我们可以查看具体的运行情况，来查找错误 - Log.final.out
它包含的是对比完以后的对比统计信息 - SJ.out.tab
它包含了剪切的信息
nohup STAR --runThreadN 20 --genomeDir /home/data/cheng/star_index --readFilesCommand gunzip -c --readFilesIn /home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.2.fq.gz --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts &
nohup STAR --runThreadN 4 --genomeDir /home/data/cheng/star_index --readFilesCommand gunzip -c --readFilesIn /home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/NALM6-ctrl_combined_R1.filtered.2.fq.gz --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix NALM6-ctrl &
nohup STAR --runThreadN 4 --genomeDir /home/data/cheng/star_index --readFilesCommand gunzip -c --readFilesIn /home/data/cheng/project/raw/NALM6-treat_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/NALM6-treat_combined_R1.filtered.2.fq.gz --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix NALM6-treat &
nohup STAR --runThreadN 4 --genomeDir /home/data/cheng/star_index --readFilesCommand gunzip -c --readFilesIn /home/data/cheng/project/raw/RS411-ctrl_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/RS411-ctrl_combined_R1.filtered.2.fq.gz --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix RS411-ctrl &
nohup STAR --runThreadN 4 --genomeDir /home/data/cheng/star_index --readFilesCommand gunzip -c --readFilesIn /home/data/cheng/project/raw/RS411-treat_combined_R1.filtered.1.fq.gz /home/data/cheng/project/raw/RS411-treat_combined_R1.filtered.2.fq.gz --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts --outFileNamePrefix RS411-treat &
# Move the BAM file into the correct folder
mv -v *Aligned.sortedByCoord.out.bam /home/data/cheng/project/results/4_aligned_sequences/aligned_bam/
# Move the logs into the correct folder
mv -v *.out /home/data/cheng/project/results/4_aligned_sequences/aligned_logs/
Step 5. Summarizing Gene Counts with featureCounts
https://www.ncbi.nlm.nih.gov/pubmed/24227677
Description
"featureCounts is a highly efficient general-purpose read summarization program that counts mapped reads for genomic features such as genes, exons, promoter, gene bodies, genomic bins and chromosomal locations. It can be used to count both RNA-seq and genomic DNA-seq reads. featureCounts takes as input SAM/BAM files and an annotation file including chromosomal coordinates of features. It outputs numbers of reads assigned to features (or meta-features). It also outputs stat info for the overall summrization results, including number of successfully assigned reads and number of reads that failed to be assigned due to various reasons (these reasons are included in the stat info)."
Now that we have our .BAM alignment files, we can then proceed to try and summarize these coordinates into genes and abundances. To do this we must summarize the reads using featureCounts or any other read summarizer tool, and produce a table of genes by samples with raw sequence abundances. This table will then be used to perform statistical analysis and find differentially expressed genes.
# Help
featureCounts -h
# Change directory into the aligned .BAM folder
cd /home/data/cheng/project/results/4_aligned_sequences/aligned_bam
# Store list of files as a variable
dirlist=$(ls -t ./*.bam | tr '\n' ' ')
echo $dirlist
# Run featureCounts on all of the samples (~10 minutes)
featureCounts \
-a /home/data/reference/gtf/* \
-o /home/data/cheng/project/results/5_final_counts/final_counts.txt \
-g 'gene_name' \
-T 4 \
$dirlist
(rna) cheng@super-Super-Server:/home/data/cheng/project/results/4_aligned_sequences/aligned_bam$ featureCounts -a /home/data/reference/gtf/* -o /home/data/cheng/project/results/5_final_counts/final_counts.txt -g 'gene_name' -T 4 $dirlist
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
v2.0.1
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 4 BAM files ||
|| o NALM6-treatAligned.sortedByCoord.out.bam ||
|| o RS411-ctrlAligned.sortedByCoord.out.bam ||
|| o RS411-treatAligned.sortedByCoord.out.bam ||
|| o NALM6-ctrlAligned.sortedByCoord.out.bam ||
|| ||
|| Output file : final_counts.txt ||
|| Summary : final_counts.txt.summary ||
|| Annotation : gencode.v35.annotation.gtf (GTF) ||
|| Dir for temp files : /home/data/cheng/project/results/5_final_counts ||
|| ||
|| Threads : 4 ||
|| Level : meta-feature level ||
|| Paired-end : no ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
\\============================================================================//
//================================= Running ==================================\\
|| ||
|| Load annotation file gencode.v35.annotation.gtf ... ||
|| Features : 1398443 ||
|| Meta-features : 59609 ||
|| Chromosomes/contigs : 25 ||
|| ||
|| Process BAM file NALM6-treatAligned.sortedByCoord.out.bam... ||
|| WARNING: Paired-end reads were found. ||
|| Total alignments : 86954860 ||
|| Successfully assigned alignments : 37685067 (43.3%) ||
|| Running time : 0.57 minutes ||
|| ||
|| Process BAM file RS411-ctrlAligned.sortedByCoord.out.bam... ||
|| WARNING: Paired-end reads were found. ||
|| Total alignments : 69733102 ||
|| Successfully assigned alignments : 30363294 (43.5%) ||
|| Running time : 0.46 minutes ||
|| ||
|| Process BAM file RS411-treatAligned.sortedByCoord.out.bam... ||
|| WARNING: Paired-end reads were found. ||
|| Total alignments : 81999795 ||
|| Successfully assigned alignments : 34018798 (41.5%) ||
|| Running time : 0.56 minutes ||
|| ||
|| Process BAM file NALM6-ctrlAligned.sortedByCoord.out.bam... ||
|| WARNING: Paired-end reads were found. ||
|| Total alignments : 70082894 ||
|| Successfully assigned alignments : 32230367 (46.0%) ||
|| Running time : 0.48 minutes ||
|| ||
|| Write the final count table. ||
|| Write the read assignment summary. ||
|| ||
|| Summary of counting results can be found in file "/home/data/cheng/projec ||
|| t/results/5_final_counts/final_counts.txt.summary" ||
|| ||
\\============================================================================//
下游分析
four_sample_deg
Cheng Liangping
24 十月, 2020
1 Importing Gene Counts into R/RStudio
Once the workflow has completed, you can now use the gene count table as an input into DESeq2 for statistical analysis using the R-programming language. It is highly reccomended to use RStudio when writing R code and generating R-related analyses. You can download RStudio for your system here: https://www.rstudio.com/products/rstudio/download/
1.1 Install required R-libraries
rm(list = ls())
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos
## CRAN
## "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"
options()$BioC_mirror
## [1] "http://mirrors.tuna.tsinghua.edu.cn/bioconductor/"
# install.packages("gage")
# BiocManager::install("org.Hs.eg.db",ask = F,update = F)
library(DESeq2)
library(ggplot2)
library(clusterProfiler)
library(biomaRt)
library(ReactomePA)
library(DOSE)
library(KEGG.db)
library(org.Hs.eg.db)
library(pheatmap)
library(genefilter)
library(RColorBrewer)
library(GO.db)
library(topGO)
library(dplyr)
library(gage)
library(ggsci)
1.1.1 Import featureCounts output
One you have an R environment appropriatley set up, you can begin to import the featureCounts table found within the 5_final_counts folder. This tutorial will use DESeq2 to normalize and perform the statistical analysis between sample groups. Be sure to know the full location of the final_counts.txt file generate from featureCounts.
Note: If you would like to use an example final_counts.txt table, look into the example/ folder.
# Import gene counts table
# - skip first row (general command info)
# - make row names the gene identifiers
countdata <- read.table("final_counts.txt", header = TRUE, skip = 1, row.names = 1)
# Remove .bam + '..' from column identifiers
colnames(countdata) <- gsub(".bam", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub("Aligned.sortedByCoord.out", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub("..", "", colnames(countdata), fixed = T)
#save(countdata,file = "raw.Rdata")
# Remove length/char columns
countdata <- countdata[ ,c(-1:-5)]
# Make sure ID's are correct
head(countdata)
## NALM6.treat RS411.ctrl RS411.treat NALM6.ctrl
## DDX11L1 26 5 20 12
## WASH7P 1119 206 202 435
## MIR6859-1 21 6 4 9
## MIR1302-2HG 2 0 3 1
## MIR1302-2 0 0 0 0
## FAM138A 0 0 0 0
1.1.2 Import metadata text file. The SampleID’s must be the first column.
# Import metadata file
# - make row names the matching sampleID's from the countdata
metadata <- read.delim("metadata.txt", row.names = 1)
# Reorder sampleID's to match featureCounts column order.
metadata <- metadata[match(colnames(countdata), metadata$sampleid), ]
# Make sure ID's are correct
head(metadata)
## Group Replicate sampleid
## NALM6.treat Hiexp Rep2 NALM6.treat
## RS411.ctrl Loexp Rep1 RS411.ctrl
## RS411.treat Loexp Rep2 RS411.treat
## NALM6.ctrl Hiexp Rep1 NALM6.ctrl
1.1.3 Make DESeq2 object from counts and metadata
# - countData : count dataframe
# - colData : sample metadata in the dataframe with row names as sampleID's
# - design : The design of the comparisons to use.
# Use (~) before the name of the column variable to compare
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
colData = metadata,
design = ~Group)
# Find differential expressed genes
ddsMat <- DESeq(ddsMat)
1.1.4 Get basic statisics about the number of significant genes
# Get results from testing with FDR adjust pvalues
results <- results(ddsMat, pAdjustMethod = "fdr", alpha = 0.05)
# Generate summary of testing.
summary(results)
##
## out of 37196 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 2345, 6.3%
## LFC < 0 (down) : 1530, 4.1%
## outliers [1] : 0, 0%
## low counts [2] : 11169, 30%
## (mean count < 4)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Check directionality of the log2 fold changes
## Log2 fold change is set as (Loexp / Hiexp)
## Postive fold changes = Increased in Loexp
## Negative fold changes = Decreased in Loexp
mcols(results, use.names = T)
## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized counts for all samples
## log2FoldChange results log2 fold change (MLE): Group Loexp vs Hiexp
## lfcSE results standard error: Group Loexp vs Hiexp
## stat results Wald statistic: Group Loexp vs Hiexp
## pvalue results Wald test p-value: Group Loexp vs Hiexp
## padj results fdr adjusted p-values
1.2 Annotate gene symbols
After alignment and summarization, we only have the annotated gene symbols. To get more information about significant genes, we can use annoated databases to convert gene symbols to full gene names and entrez ID’s for further analysis.
1.2.1 Gather gene annotation information
# Mouse genome database (Select the correct one)
library(org.Hs.eg.db)
# Add gene full name
results$description <- mapIds(x = org.Hs.eg.db,
keys = row.names(results),
column = "GENENAME",
keytype = "SYMBOL",
multiVals = "first")
# Add gene symbol
results$symbol <- row.names(results)
# Add ENTREZ ID
results$entrez <- mapIds(x = org.Hs.eg.db,
keys = row.names(results),
column = "ENTREZID",
keytype = "SYMBOL",
multiVals = "first")
# Add ENSEMBL
results$ensembl <- mapIds(x = org.Hs.eg.db,
keys = row.names(results),
column = "ENSEMBL",
keytype = "SYMBOL",
multiVals = "first")
# Subset for only significant genes (q < 0.05)
results_sig <- subset(results, padj < 0.05)
head(results_sig)
## log2 fold change (MLE): Group Loexp vs Hiexp
## Wald test p-value: Group Loexp vs Hiexp
## DataFrame with 6 rows and 10 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## MTND1P23 2424.6057 -4.11145 0.678250 -6.06185 1.34561e-09 4.39978e-08
## MTND2P28 1813.7288 3.75228 0.738809 5.07882 3.79791e-07 8.29959e-06
## MTCO1P12 6187.1454 -3.29924 0.634548 -5.19936 1.99981e-07 4.56171e-06
## MTCO2P12 169.0384 -2.08444 0.773877 -2.69351 7.07044e-03 4.79849e-02
## MXRA8 725.5861 2.58245 0.664147 3.88837 1.00920e-04 1.26220e-03
## LINC01770 54.9444 9.28692 1.703909 5.45036 5.02681e-08 1.27891e-06
## description symbol entrez
## <character> <character> <character>
## MTND1P23 MT-ND1 pseudogene 23 MTND1P23 100887749
## MTND2P28 MT-ND2 pseudogene 28 MTND2P28 100652939
## MTCO1P12 MT-CO1 pseudogene 12 MTCO1P12 107075141
## MTCO2P12 MT-CO2 pseudogene 12 MTCO2P12 107075310
## MXRA8 matrix remodeling associated 8 MXRA8 54587
## LINC01770 long intergenic non-protein coding RNA 1770 LINC01770 102724312
## ensembl
## <character>
## MTND1P23 NA
## MTND2P28 NA
## MTCO1P12 NA
## MTCO2P12 NA
## MXRA8 ENSG00000162576
## LINC01770 NA
1.2.2 Write all the important results to .txt files
# Write normalized gene counts to a .txt file
write.table(x = as.data.frame(counts(ddsMat), normalized = T),
file = 'normalized_counts.txt',
sep = '\t',
quote = F,
col.names = NA)
# Write significant normalized gene counts to a .txt file
write.table(x = counts(ddsMat[row.names(results_sig)], normalized = T),
file = 'normalized_counts_significant.txt',
sep = '\t',
quote = F,
col.names = NA)
# Write the annotated results table to a .txt file
write.table(x = as.data.frame(results),
file = "results_gene_annotated.txt",
sep = '\t',
quote = F,
col.names = NA)
# Write significant annotated results table to a .txt file
write.table(x = as.data.frame(results_sig),
file = "results_gene_annotated_significant.txt",
sep = '\t',
quote = F,
col.names = NA)
1.3 Plotting Gene Expression Data
There are multiple ways to plot gene expression data. Below we are only listing a few popular methods, but there are many more resources (Going Further) that will walk through different R commands/packages for plotting.
1.3.1 PCA plot
# Convert all samples to rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)
# Plot PCA by column variable
plotPCA(ddsMat_rlog, intgroup = "Group", ntop = 500) +
theme_bw() + # remove default ggplot2 theme
geom_point(size = 5) + # Increase point size
scale_y_continuous(limits = c(-5, 5)) + # change limits to fix figure dimensions
ggtitle(label = "Principal Component Analysis (PCA)",
subtitle = "Top 500 most variable genes")
image.png
1.3.2 Heatmap
# Convert all samples to rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)
# Gather 30 significant genes and make matrix
mat <- assay(ddsMat_rlog[row.names(results_sig)])[1:40, ]
# Choose which column variables you want to annotate the columns by.
annotation_col = data.frame(
Group = factor(colData(ddsMat_rlog)$Group),
Replicate = factor(colData(ddsMat_rlog)$Replicate),
row.names = colData(ddsMat_rlog)$sampleid
)
# Specify colors you want to annotate the columns by.
ann_colors = list(
Group = c(Loexp = "lightblue", Hiexp = "darkorange"),
Replicate = c(Rep1 = "darkred", Rep2 = "forestgreen")
)
# Make Heatmap with pheatmap function.
## See more in documentation for customization
pheatmap(mat = mat,
color = colorRampPalette(brewer.pal(9, "YlOrBr"))(255),
scale = "row", # Scale genes to Z-score (how many standard deviations)
annotation_col = annotation_col, # Add multiple annotations to the samples
annotation_colors = ann_colors,# Change the default colors of the annotations
fontsize = 6.5, # Make fonts smaller
cellwidth = 55, # Make the cells wider
show_colnames = F)
image.png
1.3.3 Volcano Plot
# Gather Log-fold change and FDR-corrected pvalues from DESeq2 results
## - Change pvalues to -log10 (1.3 = 0.05)
data <- data.frame(gene = row.names(results),
pval = -log10(results$padj),
lfc = results$log2FoldChange)
# Remove any rows that have NA as an entry
data <- na.omit(data)
# Color the points which are up or down
## If fold-change > 0 and pvalue > 1.3 (Increased significant)
## If fold-change < 0 and pvalue > 1.3 (Decreased significant)
data <- mutate(data, color = case_when(data$lfc > 0 & data$pval > 1.3 ~ "Increased",
data$lfc < 0 & data$pval > 1.3 ~ "Decreased",
data$pval < 1.3 ~ "nonsignificant"))
# Make a basic ggplot2 object with x-y values
vol <- ggplot(data, aes(x = lfc, y = pval, color = color))
# Add ggplot2 layers
vol +
ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") +
geom_point(size = 2.5, alpha = 0.8, na.rm = T) +
scale_color_manual(name = "Directionality",
values = c(Increased = "#008B00", Decreased = "#CD4F39", nonsignificant = "darkgray")) +
theme_bw(base_size = 14) + # change overall theme
theme(legend.position = "right") + # change the legend
xlab(expression(log[2]("Loexp" / "Hiexp"))) + # Change X-Axis label
ylab(expression(-log[10]("adjusted p-value"))) + # Change Y-Axis label
geom_hline(yintercept = 1.3, colour = "darkgrey") + # Add p-adj value cutoff line
scale_y_continuous(trans = "log1p") # Scale yaxis due to large p-values
image.png
1.3.4 MA Plot
https://en.wikipedia.org/wiki/MA_plot
plotMA(results, ylim = c(-5, 5))
image.png
1.3.5 Plot Dispersions
plotDispEsts(ddsMat)
image.png
1.3.6 Single gene plot
# Convert all samples to rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)
# Get gene with highest expression
top_gene <- rownames(results)[which.min(results$log2FoldChange)]
# Plot single gene
plotCounts(dds = ddsMat,
gene = top_gene,
intgroup = "Group",
normalized = T,
transform = T)
image.png
1.4 Finding Pathways from Differential Expressed Genes
Pathway enrichment analysis is a great way to generate overall conclusions based on the individual gene changes. Sometimes individiual gene changes are overwheling and are difficult to interpret. But by analyzing the pathways the genes fall into, we can gather a top level view of gene responses. You can find more information about clusterProfiler here: http://bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html
1.4.1 Set up matrix to take into account EntrezID’s and fold changes for each gene
# Remove any genes that do not have any entrez identifiers
results_sig_entrez <- subset(results_sig, is.na(entrez) == FALSE)
# Create a matrix of gene log2 fold changes
gene_matrix <- results_sig_entrez$log2FoldChange
# Add the entrezID's as names for each logFC entry
names(gene_matrix) <- results_sig_entrez$entrez
# View the format of the gene matrix
##- Names = ENTREZ ID
##- Values = Log2 Fold changes
head(gene_matrix)
## 100887749 100652939 107075141 107075310 54587 102724312
## -4.111451 3.752277 -3.299241 -2.084445 2.582450 9.286916
gene_matrix=sort(gene_matrix,decreasing = T)
1.4.2 Enrich genes using the KEGG database
kegg_enrich <- enrichKEGG(gene = names(gene_matrix),
organism = 'human',
pvalueCutoff = 0.05,
qvalueCutoff = 0.10)
kegg_enrich <- setReadable(kegg_enrich, org.Hs.eg.db, keyType = "ENTREZID")
# Plot results
barplot(kegg_enrich,
drop = TRUE,
showCategory = 10,
title = "KEGG Enrichment Pathways",
font.size = 8)
image.png
1.4.3 Enrich genes using the Gene Onotlogy
go_enrich <- enrichGO(gene = names(gene_matrix),
OrgDb = 'org.Hs.eg.db',
readable = T,
ont = "BP",
pvalueCutoff = 0.05,
qvalueCutoff = 0.10)
# Plot results
barplot(go_enrich,
drop = TRUE,
showCategory = 10,
title = "GO Biological Pathways",
font.size = 8)
image.png
1.5 Plotting KEGG Pathways
Pathview is a package that can take KEGG identifier and overlay fold changes to the genes which are found to be significantly different. Pathview also works with other organisms found in the KEGG database and can plot any of the KEGG pathways for the particular organism.
# Load pathview
library(pathview)
## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
##
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################
# Plot specific KEGG pathways (with fold change)
## pathway.id : KEGG pathway identifier
pathview(gene.data = gene_matrix,
pathway.id = "00020",
species = "hsa")
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/cheng/Downloads/TYL1/20201023
## Info: Writing image file hsa00020.pathview.png
1.6 疾病通路
有一显著差异基因富集到一个已经疾病通路： 相关疾病为：Lewy body dementia。
x <- enrichDO(gene = names(gene_matrix),
ont = "DO",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
minGSSize = 5,
maxGSSize = 500,
qvalueCutoff = 0.05,
readable = FALSE)
x <- setReadable(x, org.Hs.eg.db, keyType = "ENTREZID")
head(x)
## ID Description GeneRatio
## DOID:0060049 DOID:0060049 autoimmune disease of urogenital tract 26/1207
## DOID:12236 DOID:12236 primary biliary cirrhosis 26/1207
## DOID:1037 DOID:1037 lymphoblastic leukemia 99/1207
## DOID:0060100 DOID:0060100 musculoskeletal system cancer 98/1207
## DOID:9119 DOID:9119 acute myeloid leukemia 33/1207
## DOID:16 DOID:16 integumentary system disease 88/1207
## BgRatio pvalue p.adjust qvalue
## DOID:0060049 73/8007 1.097364e-05 0.005426588 0.004509630
## DOID:12236 73/8007 1.097364e-05 0.005426588 0.004509630
## DOID:1037 443/8007 1.839793e-05 0.005426588 0.004509630
## DOID:0060100 439/8007 2.127302e-05 0.005426588 0.004509630
## DOID:9119 107/8007 2.550088e-05 0.005426588 0.004509630
## DOID:16 394/8007 5.539968e-05 0.009824211 0.008164164
## geneID
## DOID:0060049 SNCA/CAV1/NOS2/ICOS/ABCB4/PTPN13/IL7R/ENTPD1/ITGAL/TGFBR2/TNFSF13/TNFSF13B/TLR4/FN1/TNFSF10/TNFSF9/TGFB1/HLA-DRB1/FAS/ICAM1/HGF/CXCL10/CD27/CNR1/FCRL3/ENTPD2
## DOID:12236 SNCA/CAV1/NOS2/ICOS/ABCB4/PTPN13/IL7R/ENTPD1/ITGAL/TGFBR2/TNFSF13/TNFSF13B/TLR4/FN1/TNFSF10/TNFSF9/TGFB1/HLA-DRB1/FAS/ICAM1/HGF/CXCL10/CD27/CNR1/FCRL3/ENTPD2
## DOID:1037 NRP1/FLT3/CD44/ERG/CRY1/LAIR1/HSPB1/ADRB2/SLIT2/BLACE/FCMR/ZAP70/FOSL2/IL24/TNFRSF11A/FMOD/NOS2/WT1/CDX2/BGN/IL15/CD34/RAG1/NT5E/F13A1/CSPG4/IGFBP2/MDK/IRF4/IL1B/PTPN6/FLT4/ANPEP/SLC29A2/SDC1/CD33/RUNX1/ZFP36L1/MEF2C/TP73/ARID5B/TNFSF13/FCGR2A/TNFSF13B/LILRB1/CD9/EPHB4/JAK3/CDK6/LMO2/LEF1/CISH/RB1/BCL2/RAPGEF3/CCND3/MYB/DMD/ZMIZ1/CDKN1B/TNFRSF13C/POU2AF1/IGHD/TNFSF10/TGFB1/STAT1/HLA-DRB1/IRF8/BACH2/FAS/CD40/PDE4B/CD83/FCRL5/HGF/LAMA5/FCER2/HLF/ASNS/PRKCB/FCRL2/TNFRSF8/IKZF3/PEG10/BCL2A1/TP63/BIRC3/HES1/TTC12/EFNA5/RHD/FCRL3/MME/CDKN2B/TNFRSF10A/AICDA/MLLT3/CCL22/CDKN2A
## DOID:0060100 CCN2/PDPN/CD44/ERG/HMGA2/HSPB1/NOG/PROM1/NDRG1/SPRY1/ANGPT1/SLIT2/TIMP1/BMP2/ADAM28/CAV1/LPAR1/F2RL3/FBN2/S100A6/TNFRSF11A/NOS2/HPSE/WT1/NES/CDH11/SPP1/CD34/FSCN1/SALL2/PDGFA/RUNX2/F2R/MCAM/PDGFRB/IL1B/HBEGF/JUP/LGALS1/PLAU/ERBB2/PLAUR/TLE1/KLRK1/CD33/PRUNE2/S100A4/CCL5/PGF/CREB3L1/TP73/TNFSF13B/CD9/CD99/FBN1/DUSP1/FN1/RB1/CAPN2/BCL2/DMD/RASSF1/CDKN1B/IGF2BP2/ASS1/FLI1/TNFRSF10B/PRKCA/URGCP/DYRK1B/OBSCN/TGFB1/LRP1/GLI1/ITGB3/ID1/FAS/ICAM1/HGF/UTS2R/SMO/PRKCG/PAX2/CDH1/TP63/ACKR3/EGFR/IGFBP6/DCN/LUM/GFAP/DNASE1L3/EXT1/CDKN2B/WNT5A/UCHL1/MTAP/CDKN2A
## DOID:9119 NRP1/MEIS1/CD96/PLA2G4A/CAV1/SALL4/LYZ/WT1/CDX2/RAG1/IGFBP2/IL1B/PLAUR/SELL/ITGAL/RUNX1/MGMT/IL3RA/LMO2/RB1/BCL2/CDKN1B/LYL1/IRF8/FAS/CBFA2T3/SKI/RUNX3/CAMP/CDKN2B/TNFRSF10A/KCNH2/CDKN2A
## DOID:16 PDPN/LAMA3/HSPB1/HR/ADRB2/EDA2R/ADAM33/CCL28/FLG/RIN2/SLC16A2/HRH4/TIMP1/VIM/CAV1/CD93/IL24/SERPING1/GHR/CCL1/NLRP3/IL15/CD34/NPY/COL6A5/IL1RL1/IL9R/RETN/F2R/EPX/IRAK3/IL1B/XPNPEP2/IL7R/TBC1D4/PLAU/PLAUR/SELL/ANPEP/SLC29A2/COL7A1/NR3C2/GPX1/NTRK1/TNFRSF1B/CCL5/PGF/TNFSF13/TNFSF13B/TLR4/NOS3/BCL2/SLC29A3/NOP53/HLA-DQA1/TGFB1/HLA-DRB1/ITGB3/ITGA6/FAS/CD40/TNFSF8/ALOX5AP/ICAM1/CXCL10/LAMA5/CYP1A2/C3/ADM2/CCL3L3/CAMP/BCL2A1/TP63/CCL17/EGFR/FBLN5/LMNA/CX3CL1/ALOX5/APP/TNFRSF14/FCRL3/ST14/MME/TRPC1/CDKN2B/CCL22/CDKN2A
## Count
## DOID:0060049 26
## DOID:12236 26
## DOID:1037 99
## DOID:0060100 98
## DOID:9119 33
## DOID:16 88
1.7 GO数据可视化：
1.7.1 barplot、dotplot
用散点图展示富集到的GO terms
library(enrichplot)
library(clusterProfiler)
library(patchwork)
#(1)条带图
barplot(go_enrich,showCategory=10)
image.png
#(2)气泡图
dotplot(go_enrich)
image.png
1.7.1.1 GO有向无环图
调用topGO来实现GO有向无环图的绘制,矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。
plotGOgraph(go_enrich)
image.png
## $dag
## A graphNEL graph with directed edges
## Number of Nodes = 48
## Number of Edges = 83
##
## $complete.dag
## [1] "A graph with 48 nodes."
1.7.1.2 emapplot、goplot、cnetplot
#Gene-Concept Network
cnetplot(go_enrich, categorySize="pvalue", foldChange=names(gene_matrix),colorEdge = TRUE)
image.png
cnetplot(go_enrich, foldChange=names(gene_matrix), circular = TRUE, colorEdge = TRUE)
image.png
#Enrichment Map
emapplot(go_enrich)
image.png
#(4)展示通路关系
goplot(go_enrich)
image.png
#(5)Heatmap-like functional classification
heatplot(go_enrich,foldChange = gene_matrix)
image.png
#太多基因就会糊。可通过调整比例或者减少基因来控制。
1.7.1.3 04.upsetplot
upsetplot(go_enrich)
image.png
1.7.2 用所有有entrez的基因做GSEA分析
# Remove any genes that do not have any entrez identifiers
results_entrez <- subset(results, is.na(entrez) == FALSE)
# Create a matrix of gene log2 fold changes
gene_matrix <- results_entrez$log2FoldChange
# Add the entrezID's as names for each logFC entry
names(gene_matrix) <- results_entrez$entrez
gene_matrix=sort(gene_matrix,decreasing = T)
ego_GSE_bp <- gseGO(
geneList = gene_matrix,
OrgDb = org.Hs.eg.db,
ont = "MF",
nPerm = 1000,
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)
1.7.3 ridgeplot
#gsea
ridgeplot(ego_GSE_bp)
image.png
1.7.4 gseaplot、gseaplot2及多个geneset验证结果
#gesa图
p15 <- gseaplot(ego_GSE_bp, geneSetID = 1, by = "runningScore",
title = ego_GSE_bp$Description[1])
p16 <- gseaplot(ego_GSE_bp, geneSetID = 1, by = "preranked",
title = ego_GSE_bp$Description[1])
p17 <- gseaplot(ego_GSE_bp, geneSetID = 1, title = ego_GSE_bp$Description[1])
#还有一种展示方式
p18=gseaplot2(ego_GSE_bp, geneSetID = 1, title = ego_GSE_bp$Description[1])
cowplot::plot_grid(p15, p16, p17,p18, ncol=2, labels=LETTERS[1:6])
image.png
#多条、加table(加p值)
gseaplot2(ego_GSE_bp, geneSetID = 1:3, pvalue_table = TRUE)
image.png
#上下3联(多个geneset验证结果)
pp <- lapply(1:3, function(i) {
anno <- ego_GSE_bp[i, c("NES", "pvalue", "p.adjust")]
lab <- paste0(names(anno), "=", round(anno, 3), collapse="\n")
gsearank(ego_GSE_bp, i, ego_GSE_bp[i, 2]) + xlab(NULL) +ylab(NULL) +
annotate("text", 0, ego_GSE_bp[i, "enrichmentScore"] * .9, label = lab, hjust=0, vjust=0)
})
plot_grid(plotlist=pp, ncol=1)
image.png
1.7.5 KEGG数据GSEA分析
kk_gse <- gseKEGG(geneList = gene_matrix,
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.9,
verbose = FALSE)
#使用了所有基因进行KEGG的GSEA分析发现有上调16个通路得到验证，下调0个通路显著。
down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可视化
source("kegg_plot_function.R")
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
image.png
1.8 查看kegg富集通路及差异基因可视化
1.8.1 Upsetplot for gene set enrichment analysis.
library(enrichplot)
library(clusterProfiler)
upsetplot(kegg_enrich)
image.png
1.8.2 可通过浏览器打开kegg网站查看具体信号通路的信息
head(kegg_enrich)[,1:6]
## ID Description GeneRatio
## hsa04015 hsa04015 Rap1 signaling pathway 56/1085
## hsa04060 hsa04060 Cytokine-cytokine receptor interaction 68/1085
## hsa04640 hsa04640 Hematopoietic cell lineage 30/1085
## hsa04390 hsa04390 Hippo signaling pathway 39/1085
## hsa05146 hsa05146 Amoebiasis 28/1085
## hsa04672 hsa04672 Intestinal immune network for IgA production 17/1085
## BgRatio pvalue p.adjust
## hsa04015 210/8076 1.830708e-07 0.0000583996
## hsa04060 295/8076 3.267347e-06 0.0005211418
## hsa04640 99/8076 8.949786e-06 0.0009516606
## hsa04390 157/8076 7.554828e-05 0.0060249750
## hsa05146 102/8076 1.251732e-04 0.0066839043
## hsa04672 49/8076 1.257161e-04 0.0066839043
#browseKEGG(kk.diff, 'hsa04015')
dotplot(kegg_enrich)
image.png
barplot(kegg_enrich,showCategory=10)
image.png
cnetplot(kegg_enrich, categorySize="pvalue", foldChange=gene_matrix,colorEdge = TRUE)
image.png
cnetplot(kegg_enrich, foldChange=gene_matrix, circular = TRUE, colorEdge = TRUE)
image.png
emapplot(kegg_enrich)
image.png
heatplot(kegg_enrich,foldChange = gene_matrix)
image.png
1.8.3 Going further with RNAseq analysis
You can the links below for a more in depth walk through of RNAseq analysis using R:
- http://www.bioconductor.org/help/workflows/rnaseqGene/
- http://bioconnector.org/workshops/r-rnaseq-airway.html
- http://www-huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html
- http://www-huber.embl.de/users/klaus/Teaching/DESeq2.pdf
- https://web.stanford.edu/class/bios221/labs/rnaseq/lab_4_rnaseq.html
- http://www.rna-seqblog.com/which-method-should-you-use-for-normalization-of-rna-seq-data/
- http://www.rna-seqblog.com/category/technology/methods/data-analysis/data-visualization/
- http://www.rna-seqblog.com/category/technology/methods/data-analysis/pathway-analysis/
- http://www.rna-seqblog.com/inferring-metabolic-pathway-activity-levels-from-rna-seq-data/
1.8.3.1 Citations:
-
Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc
-
Martin, Marcel. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, [S.l.], v. 17, n. 1, p. pp. 10-12, may. 2011. ISSN 2226-6089. Available at: http://journal.embnet.org/index.php/embnetjournal/article/view/200. doi:http://dx.doi.org/10.14806/ej.17.1.200.
-
Kopylova E., Noé L. and Touzet H., “SortMeRNA: Fast and accurate filtering of ribosomal RNAs in metatranscriptomic data”, Bioinformatics (2012), doi: 10.1093/bioinformatics/bts611
-
Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21. doi:10.1093/bioinformatics/bts635.
-
Lassmann et al. (2010) “SAMStat: monitoring biases in next generation sequencing data.” Bioinformatics doi:10.1093/bioinformatics/btq614 [PMID: 21088025]
-
Liao Y, Smyth GK and Shi W (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7):923-30.
-
Love MI, Huber W and Anders S (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology, 15, pp. 550.
-
Yu G, Wang L, Han Y and He Q (2012). “clusterProfiler: an R package for comparing biological themes among gene clusters.” OMICS: A Journal of Integrative Biology, 16(5), pp. 284-287.
-
Philip Ewels, Måns Magnusson, Sverker Lundin and Max Käller. “MultiQC: Summarize analysis results for multiple tools and samples in a single report” Bioinformatics (2016). doi: 10.1093/bioinformatics/btw354. PMID: 27312411
