学徒第2周是RNA-seq数据分析实战训练，讲义大纲文末的阅读原文，配套视频在B站：
九月学徒已经结业，表现还不错，学了几个NGS组学数据处理加上部分单细胞，随机安排的文献数据处理图表复现也完成的还不赖，昨天在生信技能树的WGCNA代码就是他写的；重复一篇WGCNA分析的文章（代码版）
说实话，我看到他上交的转录组学习成果其实很头疼，几十张图片排版到微信公众号真的是力气活，上一次做这件事还是：七步走纯R代码通过数据挖掘复现一篇实验文章（第1到6步）排版花了半个小时！
所以这次是我逼着学徒自己排版的，反正这辈子就苦这么一回！“作词作曲”都是他自己！
一、原始数据的下载
1.先下载所有样本的SRR号的文件
下载后得到一个SRR_Acc_List.txt文件。里面有该实验的每个样本的SRR号。
将文件上传到服务器上。放到/project/home/lyang/sra/GSE130398/1.sra_data下。
2.使用sratoolkit 的prefetch 功能下载SRA 数据
##在sra路径下新建一个文件夹1.sra_data用来存放即将下载的sra数据。
mkdir -p /project/home/lyang/sra/GSE130398/1.sra_data
##依次将SRR_Acc_List.txt中的SRR号赋值给变量id###-O 设置输出目录，默认是当前文件夹
cat SRR_Acc_List_GSE130398.txt | while read id
do
echo prefetch ${id} -O /project/home/lyang/sra/GSE130398/1.sra_data/
done > prefetch.command
##放到后台去下载
nohup bash prefetch.command &
PS: 这样的把每个样本的命令存放在 prefetch.command 脚本里面并不是我教的！
3.原始SRA文件转格式为fq文件
该过程比较耗费时间。无法设置线程数来加速转换。下次可以考虑同时并行多个xshell窗口来同时处理文件。
首先需要找到文库是双端还是单端测序：显示为双端测序
##做一个软连接文件
mkdir /project/home/lyang/sra/GSE130398/2.raw_fq
ln -s /project/home/lyang/sra/GSE130398/1.sra_data/* /project/home/lyang/sra/GSE130398/2.raw_fq/
cd /project/home/lyang/sra/GSE130398/2.raw_fq/
##创建文件转换fastq.command脚本文件
for i in `ls *.sra`
do
echo "fastq-dump --gzip --split-3 -O /project/home/lyang/sra/GSE130398/2.raw_fq/ $i"
done > fastq.command
##运行脚本文件进行批量转换文件格式
nohup bash fastq.command &
SRR8980083_1.fastq.gz是一个双端测序文件，经过fastq-dump转换后形成两个文件，分别为：
SRR8980083_1.fastq.gz
SRR8980083_2.fastq.gz
PS: 做软连接挺好的，这里如果要多个样本并行，并不需要开多个xshell窗口。可以使用控制脚本，控制代码大概如下：
mkdir -p raw_fq
conda activate qc
dump=fastq-dump
cat $config_file |while read id
do
echo $id
arr=($id)
srr=${arr[1]}
sample=${arr[0]}
if((i%number2))
then
if [ ! -f ok.dump.$srr.status ]; then
sample -O $analysis_dir --gzip --split-3 $srr.sra
touch ok.dump.$srr.status
fi
fi
i=$((i+1))
done
假设我们有100个样本，就可以使用下面的脚本控制成为6批运行，相当于每次批量处理6个样本！
step2: convert sra files to fastq files
for i in {0..5};do bash step0-sra2fastq.sh raw_fq config.sra 6 $i;done
这里其实是提交了6个脚本！
不过，一般来说，大家的服务器是有任务调度系统的，很有可能是用不上这个脚本，我这里给学徒的是小型服务器，并没有安装复杂的任务调度系统。
二、测序数据的质量控制
1.对测序结果进行测序质量统计
##做一个软连接文件
mkdir /project/home/lyang/sra/GSE130398/3.fastq_qc
ln -s /project/home/lyang/sra/GSE130398/2.raw_fq/*.fastq.gz /project/home/lyang/sra/GSE130398/3.fastq_qc/
##生成的fastqc放到~/sra/3.fastq_qc/中,-t指定线程数。
fastqc -t 10 -o /project/home/lyang/sra/GSE130398/3.fastq_qc /project/home/lyang/sra/GSE130398/3.fastq_qc/*.fastq.gz
##使用MutliQC整合FastQC结果。###注意这里是将后缀为.zip的文件进行multiqc处理
multiqc /project/home/lyang/sra/GSE130398/3.fastq_qc/*zip -o /project/home/lyang/sra/GSE130398/3.fastq_qc/
关于MiltiQC报告：
整合了所以文件的Fastqc报告，查看起来非常方便。
可以看到部分序列的接头还存在。
对其中一个文件进行查看，得知其使用的是Illumina 1.9，这对后续trim_galore操作有指导意义。
2.使用trim_galore进行质量控制
注意：用conda安装trim_galore时，名称写为trim-galore，在使用时写为trim_galore
##做一个软连接文件,将格式转换后的fq.gz文件链接至此
mkdir /project/home/lyang/sra/GSE130398/4.trim_galore
cd /project/home/lyang/sra/GSE130398/4.trim_galore
ln -s /project/home/lyang/sra/GSE130398/2.raw_fq/*.fastq.gz /project/home/lyang/sra/GSE130398/4.trim_galore/
#########i=${i/_1.fastq.gz/}意思是除去“i”后面的“_1.fastq.gz”
for i in `ls *_1.fastq.gz`
do
i=${i/_1.fastq.gz/}
echo "trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired -o /project/home/lyang/sra/GSE130398/4.trim_galore ${i}_1.fastq.gz ${i}_2.fastq.gz"
done > trim_galore.command
echo "multiqc /project/home/lyang/sra/GSE130398/4.trim_galore/*zip -o /project/home/lyang/sra/GSE130398/4.trim_galore/" >> trim_galore.command
bash trim_galore.command
做完修剪后的文件：
原文件
SRR8980083_1.fastq.gz
SRR8980083_2.fastq.gz
修剪后生成的文件
SRR8980083_1_val_1.fq.gz###修剪完成后的目标文件
SRR8980083_1_trimmed.fq.gz###中间文件，最后会被删除
SRR8980083_1.fastq.gz_trimming_report.txt###中间文件，最后会被删除
最后的文件
SRR8980083_2_val_2.fq.gz
SRR8980083_2_trimmed.fq.gz
SRR8980083_2.fastq.gz_trimming_report.txt
修剪后的fastQC质控情况如下：
接头序列被完全清除干净。
三、使用hisat2将转录组数据比对到参考基因组
1.索引的构建
在mapping过程中，使用的参考基因组都是同一个文件，如hg38.fa，如果要比对到基因组和比对到转录组，只是选择使用的软件不同，不同的软件可以对参考序列有不同的处理。对于比对到转录组的软件，它可以自己分析并去除内含子之类的序列，然后再进行比对。
方法一：官网构建好的索引地址：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/，官网上有的最好下载使用。
方法二：自行构建索引：自行下载hg38.fa文件，使用hisat2-build命令构建基因组索引
我自己下载了hg38.fa进行索引的构建。
hisat2-build hg38.fa genome
2.使用HISAT2对转录组数据进行比对
目前mapping的工具有很多，比如bwa, hisat, star等。hisat2 是其中速度最快的。同时支持DNA和RNA数据的比对。
hisat2输出的比对好的sam文件，可以通过管道无缝连接转为bam格式，以及排序，也可以分开进行。
##做一个软连接文件,将修剪后的file_trimmed.fq.gz文件链接至此，现在链接的不再是从sra转换来的fq文件了。
##之前从sra转换来的fq文件后缀是.fastq.gz，现在用的后缀为.fq.gz
mkdir /project/home/lyang/sra/GSE130398/5.mapping
cd /project/home/lyang/sra/GSE130398/5.mapping
ln -s /project/home/lyang/sra/GSE130398/4.trim_galore/*.fq.gz /project/home/lyang/sra/GSE130398/5.mapping/
##hisat2比对，生成sam文件，并且转换为bam文件。
index=/project/home/lyang/refdata/hisat/human/hg38/genome
for i in `ls *_1_val_1.fq.gz`
do
i=${i/_1_val_1.fq.gz/}
echo "hisat2 -p 10 -x ${index} -1 ${i}_1_val_1.fq.gz -2 ${i}_2_val_2.fq.gz -S ${i}.sam && samtools view -bS ${i}.sam > /project/home/lyang/sra/GSE130398/5.mapping/${i}.bam && samtools sort -@ 10 -o /project/home/lyang/sra/GSE130398/5.mapping/${i}.sort.bam ${i}.bam"
done > hisat2.command
##运行脚本
bash hisat2.command
## 为IGV可视化构建索引，软件也需要bam文件的索引，构建成功后将生成后缀为.bai文件。
for i in `ls *.sort.bam`
do
i=${i/.sort.bam/}
echo "samtools index ${i}.sort.bam"
done > index.command
##运行脚本
bash index.command
IGV安装时Windows电脑最好直接安装在默认目录下，我修改了目录到其他盘中后，连续尝试了2次安装，安装后均无法打开IGV的主页面。后来将IGV和Java都安装到了C盘后成功启动IGV主页面。
四、对基因表达进行定量
常用的基于比对的基因定量软件：Htseq-count，bedtools mutilcov，featureCount。
1.使用featureCount进行alignment-based的定量
featureCount是subread套件的一个模块，最大的优点就是速度非常快，使用全部overlap的reads计数，灵活考虑多比对的reads的计数。
所以在安装时应：
conda install -y subread
关于使用：
## 使用sort后的bam文件进行操作
mkdir /project/home/lyang/sra/GSE130398/6.featureCounts
cd /project/home/lyang/sra/GSE130398/6.featureCounts
ln -s project/home/lyang/sra/GSE130398/5.mapping/*.sort.bam
## reads计数，其中的-t,-g都是默认的。其中-g可以指定显示gtf文件中attributes那一列中的任意值。（例：gene_id，gene_name等）
featureCounts -T 10 -p -t exon -g gene_id -a
/project/home/lyang/refdata/gtf/human/gencode.v31.annotation.gtf.gz -o /project/home/lyang/sra/GSE130398/6.featureCounts/all.id.txt *.sort.bam
## 简化结果，去除特征列，只保留基因计数的列
cat all.id.txt | cut -f1,7- > counts.txt
2.featureCounts的结果解析
all.id.txt文件的表达矩阵：
Gene id：基因的ensemble基因号；从左到右依次：
Chr：多个外显子所在的染色体编号;
Start：多个外显子起始位点,与前面一一对应
End：多个外显子终止位点,与前面一一对应
Strand：正负链
Length：基因长度
sampleID：一列代表一个样本，数值表示比对到该基因上的read数目
3.使用salmon进行alignment-free的定量
Salmon可以快速从fastq快速得到基因表达 ，需要下载cDNA参考基因组。
构建cDNA序列的索引：下载Homo_sapiens.GRCh38.cdna.all.fa.gz 这个文件
具体代码：
##建立路径及索引
mkdir /project/home/lyang/sra/GSE130398/7.salmon
cd /project/home/lyang/sra/GSE130398/7.salmon
##http://asia.ensembl.org/info/data/ftp/index.html##找到所需物种的cdna序列链接
cp /project/home/lyang/refdata/salmaon/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz /project/home/lyang/sra/GSE130398/7.salmon
time salmon index -t Homo_sapiens.GRCh38.cdna.all.fa.gz -i hg38_salmon
约8min。
建立结果存储路径
ln *2_val_2.fq.gz *1_val_1.fq.gz /project/home/lyang/sra/GSE130398/7.salmon/
for i in `ls *_1_val_1.fq.gz`
do
i=${i/_1_val_1.fq.gz/}
echo "salmon quant -i hg38_salmon -l A -1 ${i}_1_val_1.fq.gz -2 ${i}_2_val_2.fq.gz -p 20 -o ${i}_quant"
done > salmon.command
less salmon.command
sh salmon.command
quant.sf结果文件：
name中的T表示转录本
Name：target transcript 名称，由输入的 transcript database (FASTA file)所提供。各列含义解析：
Length：target transcript 长度，即有多少个核苷酸。
EffectiveLength：target transcript 计算的有效长度。此项考虑了所有建模的因素，这将影响从这个转录本中
取样片段的概率，包括片段长度分布和序列特异性和gc片段偏好
TPM：估计转录本的表达量
NumReads：估计比对到每个转录本的reads数。
Salmon输出其他文件：
cmd_info.json：JSON格式文件，记录salmon程序运行的命令和参数
lib_format_counts.json：Observed library format counts。当运行salmon是 mapping-based mode时，则会生成改文件。JSON格式文件，记录有关文库格式和reads比对的情况。
eq_classes.txt：Equivalence class file。当Salmon运行时，应用参数--dumpEq，则会生成此文件。
aux_info：辅助文件夹，内含多个文件
fld.gz：在辅助文件夹中，该文件记录的是观察到的片段长度分布的近似值
observed_bias_3p.gz：Sequence-specific bias files
expected_gc.gz, observed_gc.gz：当Salmon运行时，应用fragment-GC bias correction，在辅助文件夹
中则会生成这两个文件。记录Fragment-GC bias。meta_info.json：JSON格式文件，记录salmon程序运行的统计信息
ambig_info.tsv：tab分隔符的文本文件，含有两列。记录的是每个转录本对应的 the number of uniquelymapping reads 和 the total number of ambiguously-mapping reads
五、转换为表达矩阵
针对featurecount的结果输出文件进行转换：
mkdir /project/home/lyang/sra/GSE130398/8.final_matrix
cp /project/home/lyang/sra/GSE130398/6.featureCounts/all.id.txt /project/home/lyang/sra/GSE130398/8.final_matrix
##去除抬头第一行，在vim里第一行按"dd"来删去第一行
vim all.id.txt
##选取有意义的矩阵信息保存到count.txt文件中。
cut -f1,7-12 --output-delimiter="," all.id.txt > count.csv
1568708315593
针对salmon的结果输出文件进行转换：
由于在R中使用函数更加灵活，故将salmon的输出文件导入R中进行处理。
由于salmon的输出文件是转录组mRNA的ID号，故需要将其转换为基因ID号。需要自行下在相关文件（基因关系文件hg38_tx2gene.txt）进行转换。
注：先在R的工作目录下新建quants文件夹，将所有的salmon输出文件复制到这个目录下。
具体R代码如下：
###使用这个脚本前，将salmon的输出文件（类似SRR8980086_quant）整个复制到R目录下的quants文件下（如果没有就新建一个）
###使用这个脚本前，将salmon的输出文件（类似SRR8980086_quant）整个复制到R目录下的quants文件下（如果没有就新建一个）
f1='hg38_tx2gene.txt'
tx2gene=read.table(f1,stringsAsFactors = F)
head(tx2gene)
library(stringr)
tx2gene[,1]=str_split(tx2gene[,1],'_',simplify = T)[,1]
tx2gene[,2]=str_split(tx2gene[,2],'_',simplify = T)[,1]
head(tx2gene)
dir=file.path(getwd(),'quants/')
dir
files <- list.files(pattern="*sf",dir,recursive=T)
files=file.path(dir,files)
all(file.exists(files))
library("tximport")
library("readr")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
names(txi)
head(txi$length)
head(txi$counts)
library(stringr)
files
sapply(strsplit(files,'\\/'), function(x) x[length(x)-1])
t1=sapply(strsplit(files,'\\/'), function(x) x[length(x)-1])
t1
gsub('_quant','',t1)
colnames(txi$counts)= gsub('_quant','',t1)
tmp=txi$counts
exprSet=apply(tmp,2,as.integer)
rownames(exprSet)=rownames(tmp)
dim(exprSet)
write.csv(exprSet,file = "salmon_exprSet.csv",row.names = T)###文件保存格式有2个，自行选择
save(exprSet,file=paste0('quants-exprSet.Rdata'))
两个表达矩阵可以简单比较一下，我们之前就写过教程：
六、R语言中对数据进行下游分析
差异分析
因为是rpkm的数据矩阵，edgeR 和 DESeq2 使用原始的count矩阵作为输入，所以这里使用limma包进行差异分析。
limma包接受多种数据类型：芯片数据，rpkm，counts（需要用voom进行normalization）。
edgeR 和 DESeq2 主要接受counts矩阵数据。
rm(list = ls())
options(stringsAsFactors = F)
rpkm <- read.table("GSE130398_fpkm.txt",header = T,sep = "\t")
colnames(rpkm)
rownames(rpkm) <- rpkm[,1]
rpkm <- rpkm[,-1]
boxplot(log(rpkm+1),outline=T)##结果看下图
group_list <- c(rep("ko",3),rep("wt",3))
exprSet <- log(rpkm+1)
g1="wt"
g2="ko"
pro='RNA_seq'
LIMMA真的作者提供的rpkm矩阵进行差异分析时候，采用了 log(rpkm+1) 形式。
library(edgeR)
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
con=paste0(g2,'-',g1)
cont.matrix=makeContrasts(contrasts=c(con),levels = design)
##step1
fit <- lmFit(exprSet, design)
##step2
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
##step3
tempOutput = topTable(fit2, coef=con, n=Inf)
DEG_limma = na.omit(tempOutput)
head(DEG_limma)
# logFC AveExpr t P.Value adj.P.Val B
# ENSG00000177302 -4.969664 7.885902 -94.54984 7.220217e-09 0.0001399139 10.304558
# ENSG00000110768 -6.631765 12.258761 -73.42904 2.358375e-08 0.0001399139 9.746754
# ENSG00000164494 -3.484015 5.780233 -72.04955 2.577333e-08 0.0001399139 9.697859
# ENSG00000089248 64.596810 81.398602 69.09276 3.135886e-08 0.0001399139 9.586232
# ENSG00000083635 -4.543421 7.143579 -68.77335 3.204644e-08 0.0001399139 9.573585
# ENSG00000257093 5.511872 7.433534 66.76276 3.682102e-08 0.0001399139 9.491156
nrDEG=DEG_limma_voom[,c(1,4)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
logFC_cutoff <- with(nrDEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)))
upgenes=rownames(DEG_limma_voom[with(DEG_limma_voom,logFC>logFC_cutoff & adj.P.Val<0.05),])
downgenes=rownames(DEG_limma_voom[with(DEG_limma_voom,logFC < -logFC_cutoff & adj.P.Val<0.05),])
##画热图及火山图
if(T){
need_DEG=nrDEG
n=paste0(pro,'_limma')
library(pheatmap)
exprSet=log10(exprSet+1)
choose_gene=head(rownames(need_DEG),100)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
g1=pheatmap(choose_matrix)
print(g1)
ggsave(g1,filename = paste0(n,'_heatmap.png'))
logFC_cutoff <- with(need_DEG,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
need_DEG$change = as.factor(ifelse(need_DEG$pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),'\nThe number of up gene is ',nrow(need_DEG[need_DEG$change =='UP',]) ,'\nThe number of down gene is ',nrow(need_DEG[need_DEG$change =='DOWN',])
)
library(ggplot2)
g2 = ggplot(data=need_DEG,aes(x=log2FoldChange, y=-log10(pvalue),color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
xlim(-5,5)+ggtitle( this_tile ) +
theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red'))
print(g2)
ggsave(g2,filename = paste0(n,'_volcano.png'))
}
save(down_gene_symbol,up_gene_symbol,file = "deg.Rdata")
参数选择标准：
找到差异基因后按padj排序，取前100个基因作图：
logFC_cutoff利用公式计算mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange))
adj.P.Val设为0.05
火山图结果：
热图结果：
KEGG的注释：
找到的差异基因，上调的779（经过去重后687 个），下调的有527（经过去重后485 个），一共有1306个（经过去重后1171个）基因，但是第一次做富集分析的时候却找不到富集的KEGG通路或者GO的结果。
代码如下：
> kk.up <- enrichKEGG(gene = gene_up,
+ organism = 'hsa',
+ #universe = gene_all,
+ pvalueCutoff = 0.05,
+ qvalueCutoff =0.2)
--> No gene can be mapped....
--> Expected input gene ID: 128,54578,7364,51703,125,2990
--> return NULL...
kk.up <- enrichKEGG(gene = gene_up,
+ organism = 'hsa',
+ #universe = gene_all,
+ pvalueCutoff = 0.9,
+ qvalueCutoff =0.9)##pvalue和qvalue设的较大，保留尽可能多的基因
--> No gene can be mapped....
--> Expected input gene ID: 5211,83401,501,84869,84532,414328
--> return NULL...
kk.down <- enrichKEGG(gene = gene_down,
+ organism = 'hsa',
+ #universe = gene_all,
+ pvalueCutoff = 0.9,
+ qvalueCutoff =0.9)
--> No gene can be mapped....
--> Expected input gene ID: 25796,5211,5105,226,223,9524
--> return NULL...
> g_list=list(gene_up=gene_up,
gene_down=gene_down,
gene_diff=gene_diff)
> go_enrich_results <- lapply( g_list , function(gene) {
lapply( c('BP','MF','CC') , function(ont) {
cat(paste('Now process ',ont ))
ego <- enrichGO(gene = gene,
#universe = gene_all,
OrgDb = org.Hs.eg.db,
ont = ont ,
pAdjustMethod = "BH",
pvalueCutoff = 0.99,
qvalueCutoff = 0.99,
readable = TRUE)
print( head(ego) )
return(ego)
})
})
Now process BP--> No gene can be mapped....
--> Expected input gene ID: 10361,54361,6406,245711,10371,9212
--> return NULL...
NULL
Now process MF--> No gene can be mapped....
--> Expected input gene ID: 80235,440275,10111,91801,79087,6419
--> return NULL...
NULL
Now process CC--> No gene can be mapped....
--> Expected input gene ID: 3066,3054,57504,1457,58516,8464
--> return NULL...
NULL
Now process BP--> No gene can be mapped....
--> Expected input gene ID: 6777,6159,3291,4436,26998,6830
--> return NULL...
NULL
Now process MF--> No gene can be mapped....
--> Expected input gene ID: 6165,26024,4361,3028,1915,29954
--> return NULL...
NULL
Now process CC--> No gene can be mapped....
--> Expected input gene ID: 23522,55869,1457,10933,10943,54815
--> return NULL...
NULL
Now process BP--> No gene can be mapped....
--> Expected input gene ID: 158880,1080,4867,3066,1312,23462
--> return NULL...
NULL
Now process MF--> No gene can be mapped....
--> Expected input gene ID: 160335,79053,5917,60678,56052,85365
--> return NULL...
NULL
Now process CC--> No gene can be mapped....
--> Expected input gene ID: 6871,5931,9898,10856,54815,6878
--> return NULL...
NULL
因为之前几乎没怎么做过KEGG或者GO，为了找原因，在网上一直查找相关教程。后来发现是因为enrichGO函数和enrichKEGG函数需要的输入id是ENTREZID，我之前用成了ENSEMBL ID。关于个各种ID之间的关系和转换，大家可以去生信技能树中查看，jimmy老师都有很详细的教程。
rm(list = ls()) ## 魔幻操作，一键清空~
load(file = 'deg.Rdata')##载入前面分析得到的差异分析结果
library(org.Hs.eg.db)
library(clusterProfiler)
tmp=toTable(org.Hs.egENSEMBL)
upgene=tmp[match(upgenes,tmp$ensembl_id),1]
downgene=tmp[match(downgenes,tmp$ensembl_id),1]
gene_up=unique(upgene)##对得到的上调基因名去重
gene_down=unique(downgene)##对得到的下调基因名去重
gene_diff=unique(c(gene_up,gene_down))####将上下调基因合并并去重
##pvalue和qvalue设的较大，保留尽可能多的基因，因为根据Y叔说可以根据得到的结果再进行自由筛选。
##具体说明可以看这里https://mp.weixin.qq.com/s/odA-xzI4lCMDmyZxtEMwFg
kk.up <- enrichKEGG(gene = upgene,
organism = 'hsa',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.up)[,1:6]
kk=kk.up
dotplot(kk)##画图初步展示
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')##直接在结果中ID注释转换
write.csv(kk@result,paste0(pro,'_kk.up.csv'))##保存本地
kk.down <- enrichKEGG(gene = downgene,
organism = 'hsa',
#universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
head(kk.down)[,1:6]
kk=kk.down
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.down.csv'))
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk.diff)[,1:6]
kk=kk.diff
dotplot(kk)
kk=DOSE::setReadable(kk, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
write.csv(kk@result,paste0(pro,'_kk.diff.csv'))
返回结果中GeneRatio一列中分母是所有注释到KEGG通路编号上的所有差异基因（并不是所有找到的差异基因都可以在通路上找到），分子是在该term中存在的差异基因数目。
BgRatio一列中分母是所有注释到KEGG通路编号上的背景基因数，分子是该term中背景基因数（不同term之间可能存在重叠）。
up上调基因
down下调基因：
diff全部差异基因：
GO注释
rm(list = ls()) ## 魔幻操作，一键清空~
load(file = 'deg.Rdata')##载入前面分析得到的差异分析结果
gene_up=unique(gene_up)##对得到的上调基因名去重
gene_down=unique(gene_down)##对得到的下调基因名去重
gene_diff=unique(c(gene_up,gene_down))####将上下调基因合并并去重
g_list=list(gene_up=gene_up,
gene_down=gene_down,
gene_diff=gene_diff)
###批量循环做go分析
go_enrich_results <- lapply( g_list , function(gene) {
lapply( c('BP','MF','CC') , function(ont) {
cat(paste('Now process ',ont ))
ego <- enrichGO(gene = gene,
#universe = gene_all,
OrgDb = org.Hs.eg.db,
ont = ont ,
pAdjustMethod = "BH",
pvalueCutoff = 0.99,
qvalueCutoff = 0.99,
readable = TRUE)
print( head(ego) )
return(ego)
})
})
##利用循环分别对上，下调基因以及所有的差异基因做'BP','MF','CC'三个方面的图，并直接保存本地
n1= c('gene_up','gene_down','gene_diff')
n2= c('BP','MF','CC')
for (i in 1:3){
for (j in 1:3){
fn=paste0(pro, '_dotplot_',n1[i],'_',n2[j],'.png')
cat(paste0(fn,'\n'))
png(fn,res=150,width = 1080)
print( dotplot(go_enrich_results[[i]][[j]] ))
dev.off()
}
}
上调差异基因：
BP
MF
CC
下调差异基因：
BP
MF
CC
上下调差异基因汇总：
BP
MF
CC
七，针对counts矩阵做差异分析
在第6步，因为作者仅仅是提供了rpkm矩阵，所以我们勉强采用了limma对 log(rpkm+1) 进行差异分析，也是得到了上下调基因及通路。不过，毕竟不是金标准，而且我们还跑了RNA-seq分析流程拿到了自己的counts矩阵，理论上我们应该是可以走真正的RNA-seq的下游分析。(还可以跟作者的rpkm矩阵比较，我们下期再见。)
一定要继续关注哦，下期更精彩！
学徒写在最后
首先感谢Jimmy老师的教程和代码，基本上只要跟着一步步学下来，肯定能复现漂亮的图，但是其中的原理需要自己仔细研究和领会。
另外，非常感谢jimmy老师对我耐心的指导和引导，当我遇到问题时，能一下了解我遇到的代码问题在哪里，比我百度谷歌半天教程都有用。
学徒已经做的很优秀了，一个月的时间总是短暂的，但学习的脚步不能停下，希望他回去以后能有更多的学习成果跟大家分享！
当然，如果你感兴趣学徒培养，也可以看看招学徒的通知：
你可以先看看我是如何培养学徒的：