九月学徒转录组学习成果展（3万字总结）（上篇） / 开普饭

学徒第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.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结果文件：

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矩阵比较，我们下期再见。)

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

学徒写在最后

首先感谢Jimmy老师的教程和代码，基本上只要跟着一步步学下来，肯定能复现漂亮的图，但是其中的原理需要自己仔细研究和领会。
另外，非常感谢jimmy老师对我耐心的指导和引导，当我遇到问题时，能一下了解我遇到的代码问题在哪里，比我百度谷歌半天教程都有用。

我写在后面

学徒已经做的很优秀了，一个月的时间总是短暂的，但学习的脚步不能停下，希望他回去以后能有更多的学习成果跟大家分享！

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