查看原文
其他


学徒第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%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报告,查看起来非常方便。

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老师对我耐心的指导和引导,当我遇到问题时,能一下了解我遇到的代码问题在哪里,比我百度谷歌半天教程都有用。

我写在后面

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


当然,如果你感兴趣学徒培养,也可以看看招学徒的通知:

你可以先看看我是如何培养学徒的:

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存