【原】九月学徒转录组学习成果展（3万字总结）（上篇）

健明 2021-07-14

展开全文

学徒第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窗口来同时处理文件。

首先需要找到文库是双端还是单端测序：显示为双端测序

1568634483954

###做一个软连接文件
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%$number1==$number2))
 then

  if [  ! -f  ok.dump.$srr.status  ]; then
   $dump -A  $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报告，查看起来非常方便。

1568888954225

1568889172606

可以看到部分序列的接头还存在。

对其中一个文件进行查看，得知其使用的是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质控情况如下：

1568889052447

1568889202244

接头序列被完全清除干净。

三、使用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

1568635031236

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主页面。

1568635104949

四、对基因表达进行定量

常用的基于比对的基因定量软件：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

1568635141624

2.featureCounts的结果解析

all.id.txt文件的表达矩阵：

1568635167599

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./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结果文件：

1568635237537

1568635283728

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'))

1568888707816

1568888719855

两个表达矩阵可以简单比较一下，我们之前就写过教程：

六、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'

1570496484215

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

火山图结果：

1570501362494

热图结果：

1570501339362

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上调基因

1570502010478

down下调基因：

1570502036244

diff全部差异基因：

1570502057961

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

1570502507741

MF

1570502468650

CC

1570502489558

下调差异基因：

BP

1570502526553

MF

1570502537347

CC

1570502550588

上下调差异基因汇总：

BP

1570502563692

MF

1570502577426

CC

1570502591835

七，针对counts矩阵做差异分析

在第6步，因为作者仅仅是提供了rpkm矩阵，所以我们勉强采用了limma对 log(rpkm+1) 进行差异分析，也是得到了上下调基因及通路。不过，毕竟不是金标准，而且我们还跑了RNA-seq分析流程拿到了自己的counts矩阵，理论上我们应该是可以走真正的RNA-seq的下游分析。(还可以跟作者的rpkm矩阵比较，我们下期再见。)

一定要继续关注哦，下期更精彩！