操作记录-2020-12-21: dongfeng_GAS_project

最新推荐文章于 2024-04-17 21:17:24 发布

垚垚爸爱学习

最新推荐文章于 2024-04-17 21:17:24 发布

阅读量950

点赞数 1

分类专栏：操作记录

本文链接：https://blog.csdn.net/xiaomotong123/article/details/111480068

版权

今天老菜鸟接了实验室一个大活，要帮学生把一个课题结掉，因此需要处理三批次数据。当然了，老菜鸟这些天正在忙一些其他事情，不知道这些事情有了结果后，还会不会再继续这些工作了，暂且记录一下自己目前复杂的心情吧。

一、RNA_seq_dongfeng_2020_12_17

1. 建立相应目录

# 新建对应的目录
mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt  scripts_log

2. 检查数据完整性

这次将提供的md5.txt中数据提前编辑了一下，修改了样品名称和md5中的目录和名称，操作简便多了。

(base) zexing@DNA:~/projects/dongfeng/RNA_seq/2020_12_17/MD5_txt$ cat md5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
../clean_data/Rescue_WT_2.clean.fq.gz: OK
../clean_data/Rescue_W93A_1.clean.fq.gz: OK
../clean_data/Rescue_WT_1.clean.fq.gz: OK
../clean_data/Rescue_W93A_2.clean.fq.gz: OK

3.在Linux服务器中对RNA_seq数据进行处理

vim新建RNA_seq_script将数据质控、比对、格式转换、排序、拼接和定量综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2020/12/21       zexing            First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/dongfeng/RNA_seq/2020_12_17

# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz

# 利用for循环进行后续操作
for i in Rescue_W93A Rescue_WT 
do
# 对数据进行比对
hisat2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_hg19/hisat2_index/hisat2_index_hg19 \
-1 ${dir}/clean_data/${i}_1.clean.fq.gz \
-2 ${dir}/clean_data/${i}_2.clean.fq.gz \
-S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 16 -G /f/xudonglab/zexing/reference/UCSC_hg19/hg19_genes.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done

后台运行RNA_seq_script：

nohup bash RNA_seq_script > RNA_seq_script_log &

4. 使用prepDE.py脚本提取read_counts数值

进入ballgown文件夹，将prepDE.py脚本拷贝至当前文件夹

cp /f/xudonglab/zexing/software/prepDE.py ./

退出当前conda环境

conda deactivate

使用python命令直接运行脚本

python prepDE.py

运行结果中"gene_count_matrix.csv"即是DESeq2的输入文件。

二、RNA_seq_dongfeng_2018_08_16

1. 建立相应目录

# 新建对应的目录
mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt  scripts_log

2. 在Linux服务器中对RNA_seq数据进行处理

vim新建RNA_seq_script将数据质控、比对、格式转换、排序、拼接和定量综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2020/12/22      zexing            First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/dongfeng/RNA_seq/2018_08_16

# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz

# 利用for循环进行后续操作
for i in KD_scr_1 KD_sh3_1 KD_sh4_1 
do
# 对数据进行比对
hisat2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_hg19/hisat2_index/hisat2_index_hg19 \
-U ${dir}/clean_data/${i}.fq.gz \
-S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 16 -G /f/xudonglab/zexing/reference/UCSC_hg19/hg19_genes.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done

后台运行RNA_seq_script：

nohup bash RNA_seq_script > RNA_seq_script_log &

3. 使用prepDE.py脚本提取read_counts数值

进入ballgown文件夹，将prepDE.py脚本拷贝至当前文件夹

cp /f/xudonglab/zexing/software/prepDE.py ./

退出当前conda环境

conda deactivate

使用python命令直接运行脚本

python prepDE.py

运行结果中"gene_count_matrix.csv"即是DESeq2的输入文件。

三、RNA_seq_dongfeng_2021_02_01

1.建立相应目录

# 新建对应的目录
mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt  scripts_log

2.检查数据完整性

cat md5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt
# 数据处理记录如下
./clean_data/W93A_DOX_1.clean.fq.gz: OK
./clean_data/WT_DOX_2.clean.fq.gz: OK
./clean_data/WT_DOX_1.clean.fq.gz: OK
./clean_data/W93A_DOX_2.clean.fq.gz: OK
md5sum: WARNING: 1 line is improperly formatted

3.在Linux服务器中对RNA_seq数据进行处理

vim新建RNA_seq_script将数据质控、比对、格式转换、排序、拼接和定量综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2021/02/01       zexing            First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/dongfeng/RNA_seq/2021_02_01

# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz

# 利用for循环进行后续操作
for i in WT_DOX W93A_DOX
do
# 对数据进行比对
hisat2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_hg19/hisat2_index/hisat2_index_hg19 \
-1 ${dir}/clean_data/${i}_1.clean.fq.gz \
-2 ${dir}/clean_data/${i}_2.clean.fq.gz \
-S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 16 -G /f/xudonglab/zexing/reference/UCSC_hg19/hg19_genes.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done

后台运行RNA_seq_script：

nohup bash RNA_seq_script > RNA_seq_script_log &

4. 使用prepDE.py脚本提取read_counts数值

进入ballgown文件夹，将prepDE.py脚本拷贝至当前文件夹

cp /f/xudonglab/zexing/software/prepDE.py ./

退出当前conda环境

conda deactivate

使用python命令直接运行脚本

python prepDE.py

运行结果中"gene_count_matrix.csv"即是DESeq2的输入文件。

四、ChIP_seq_dongfeng_2019_10_23

1.建立相应目录

# 新建对应的目录
mkdir raw_data clean_data bam bam_bw bam_sort sam macs2_bdgdiff macs2_callpeak matrix_reference_point matrix_scale_regions fastqc_report MD5_txt scripts_log

2.检查数据完整性

cat md5.txt > check_md5sum.txt && md5sum -c check_md5sum.txt

3. 在Linux服务器中对ChIP_seq数据进行处理并常规call peak。

vim新建ChIP_seq_script_1将数据质控、比对、格式转换、排序、生成目录、bamCoverage命令转换文件格式和常规callpeak综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for ChIP-seq data analysis.
# History
#     2020/12/22       zexing            First release
# 设置变量${dir}为常用目录
dir=/f/xudonglab/zexing/projects/dongfeng/ChIP_seq/2019_10_23 # 用户名称和日期需要更改
input=/f/xudonglab/zexing/reference/ChIP_input/
# 对数据进行质控
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fq.gz

# 利用for循环进行后续操作
for i in H3K18
do
# 对数据进行比对
bowtie2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_hg19/bowtie2_index/hg19 -1 ${dir}/clean_data/${i}_R1.fq.gz -2 ${dir}/clean_data/${i}_R2.fq.gz -S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据生成目录
samtools index -@ 16 ${dir}/bam_sort/${i}.bam.sort 

# bamCoverage命令转换文件格式
bamCoverage -p 16 -v -b ${dir}/bam_sort/${i}.bam.sort -o ${dir}/bam_bw/${i}.bam.sort.bw

# 使用macs2进行常规callpeak
macs2 callpeak -t ${dir}/bam_sort/${i}.bam.sort -c ${input}/U87_Normoxia_input.bam.sort -f BAM -g mm -B -q 0.05 --outdir ${dir}/macs2_callpeak/ -n ${i}

done

在后台运行ChIP_seq_script_1：

nohup bash ChIP_seq_script_1 > ChIP_seq_script_1_log &

五、作图分析

1. 对RNA_seq_2018_08_16利用DESeq2包进行差异基因分析

#This script is used for analysis of dongfeng RNA-seq data
#History
# Lizexing           2020-12-23             First release
#genecount文件来源于Stringtie软件分析，后面为本地电脑操作
# DESeq2进行差异分析 ----------------------------------------------------------------

#清空环境变量
rm(list=ls())
#设置工作目录
setwd("G:/dongfeng/RNA-seq/2018_08_16/gene_count/")
#读入基因表达值，设定行名为gene_id
gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F)
#对gene_id一列进行拆分，去除重复名称
library(stringr)
#设置空的"gene_count_1"向量，行数与上面的测序结果一致
gene_count_1<-rep(NA,nrow(gene_count))
#利用for循环，对gene_count数据框中的重复列进行拆分提取
for (i in 1:nrow(gene_count)){
   
  gene_count_1[i] <- unlist(str_split(gene_count[i,1], pattern = "\\|"))[1]
}
#显示拆分后的结果
head(gene_count_1)
#对原数据框中的特定序列重新赋值
gene_count$gene_id <- gene_count_1
#显示文件的前6行信息
head(gene_count)
#将第一列作为文件的行名
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
#显示文件的前6行信息
head(gene_count)

#将该文件保存至对应目录
write.csv(gene_count, file = "G:/dongfeng/RNA-seq/2018_08_16/gene_count/gene_count.csv", row.names = TRUE)

#加载DESeq2包
library(DESeq2)
#DESeq2需要三种矩阵，分别为countData表达矩阵,colData样品信息矩阵及design差异表达矩阵
#countData为表达矩阵即gene_count
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#设置condition样品组别、重复数
condition <- factor(c(rep("Scr", 1), rep("shGAS41", 2)), levels = c("Scr","shGAS41"))

#显示condition设置信息
condition

#设置group组对应的样品信息矩阵colData
colData <- data.frame(row.names = colnames(gene_count), condition)

#显示colData设置信息
colData

#在R里面用于构建公式对象，~左边为因变量，右边为自变量。
#标准流程：dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) 
#countData为表达矩阵即countdata
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#对dds_group进行标准流程构建
dds <- DESeqDataSetFromMatrix(gene_count, colData, design = ~condition)

#对原始dds_group进行normalize
dds <- DESeq(dds)

#显示dds信息
dds

# 对差异分析结果进行保存 -------------------------------------------------------------

#使用DESeq2包中的results()函数，提取差异分析的结果
#Usage:results(object, contrast, name, .....）
#将提取的差异分析结果定义为变量"res" 
#contrast: 定义谁和谁比较，处理组在前，对照组在后
#将group组提取分析结果并保存为res
res = results(dds, contrast=c("condition","shGAS41","Scr"))

#对结果res利用order()函数按pvalue值进行排序
#创建矩阵时，X[i,]指矩阵X中的第i行，X[,j]指矩阵X中的第j列
#order()函数先对数值排序，然后返回排序后各数值的索引，常用用法：V[order(V)]或者df[order(df$variable),]
#对res_group组进行排序
res = res[order(res$pvalue),]

#显示res结果前6行信息
head(res)

#对res_group矩阵进行总结，利用summary命令统计显示一共多少个genes上调和下调
summary(res)

#将差异分析的所有结果进行输出保存
write.csv(res, file="G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/all_different_genes/all_different_genes_genecount.csv")

#利用table函数统计显著差异基因的数目
#显著差异的定义为padj<0.05
table(res$padj<0.05)

#对具有显著性差异的结果进行过滤、提取
#获取padj小于0.05，表达倍数取以2为对数后大于1或者小于-1的差异表达基因（即表达倍数相差2倍及以上）
#使用subset()函数过滤需要的结果至新的变量significant_different_genes_group中
#Usage:subset(x, ...)，其中x为objects，...为筛选参数或条件
#对group中数据进行过滤、提取
significant_padj_different_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

#使用dim函数查看该结果的维度、规模
dim(significant_padj_different_genes)

#显示结果的前6行信息
head(significant_padj_different_genes)

#对显著差异基因进行输出保存
write.csv(significant_padj_different_genes, file = "G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/significant_padj_different_genes_genecount.csv")

注意：输出文件保留了pvalue <0.05和padj <0.05两种文件，测序中没有重复，故将两对shRNA的结果作为重复进行分析。

2. 对RNA_seq_2018_08_16绘制火山图

# 火山图 ---------------------------------------------------------------------

#代码参考网站：https://www.jianshu.com/p/e651a182c65d
#火山图的图形非常像火山喷发的形状。
#火山图通常用来展示差异表达的基因，常常出现在芯片、转录组、蛋白组、代谢组等组学检测技术的结果中，并且通常伴随热图一起出现。

#清空环境变量
rm(list=ls())
#获取当前工作目录
getwd()
#设置工作目录
setwd("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/all_different_genes")
#读取数据至deg.data变量中
#此处需要读取DESeq2分析的全部差异基因，包括显著和非显著基因
#对group_1的全部差异基因进行读取
deg.data <- read.csv("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/all_different_genes/all_different_genes_genecount.csv", header = T, sep = ",")

#更改文件行名称为第一列
rownames(deg.data) <- deg.data$X

#更改文件列名称为需要的名称
colnames(deg.data)  <- c("gene_symbol",colnames(deg.data)[c(2:7)])

#显示文件前6行查看文件信息
head(deg.data)

###画火山图只需要其中的log2FC和padj就可以，dongfeng选用pvalue值进行筛选，故此后使用pvalue进行设置
#adj.p.value为校正后的P值，因为基因和基因并不是相互独立的，所以我们需要对P值进行校正来降低结果的假阳性，常用的校正方法为FDR校正
#绘图之前需要对pvalue进行转换，可以拉开差异表达基因之间的间距
#对差异基因pvalue值进行log10转换
deg.data$logP <- -log10(deg.data$pvalue)

#开始绘制基本热图
#利用ggplot2的两个包绘制火山图
#安装ggpubr包、ggthemes包
#install.packages("ggpubr")
#install.packages("ggthemes")
#设置工作目录
setwd("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/火山图")

#加载ggpubr包
library(ggpubr)
library(ggthemes)
#x轴为实验组基因表达量比对照组基因表达量的倍数差异
#y轴则为实验组比对照组之后的pvalue值
#火山图上一个点代表一个基因，而颜色则代表他们是显著上调还是显著下调
ggscatter(deg.data, x="log2FoldChange", y="logP") + theme_base()

#上述命令出来的图很丑，需要对log2FoldChange和pvalue数据进行过滤
#新加一列Group
deg.data$Group = "not-significant"

#将pvalue<0.05且log2FC>=1的基因设为显著上调基因
deg.data$Group[which((deg.data$pvalue<0.05) & (deg.data$log2FoldChange >=1))] ="up-regulated"

#将pvalue<0.05且log2FC=<-1的基因设为显著下调基因
deg.data$Group[which((deg.data$pvalue<0.05) & (deg.data$log2FoldChange <= -1))] ="down-regulated"

#查看上调和下调基因数目
table(deg.data$Group)

#使用添加了上调和下调基因的数据重新绘制火山图
#使用color参数指定点的颜色
ggscatter(deg.data, x = "log2FoldChange", y = "logP", 
          color = "Group") + theme_base()

#修改点的大小size和更改差异表达基因的颜色palette
ggscatter(deg.data, x = "log2FoldChange", y = "logP", 
          color = "Group",
          palette = c("green", "gray", "red"),
          size = 1) + theme_base()

#使用geom_hline和geom_vline分别添加横向和纵向的辅助线
#为火山图添加logP分界线（geom_hline）和logFC分界线（geom_vline）
ggscatter(deg.data, x = "log2FoldChange", y = "logP", 
          color = "Group",
          palette = c("green", "gray", "red"),
          size = 1) + theme_base() +
  geom_hline(yintercept = 1.30, linetype="dashed") +
  geom_vline(xintercept = c(-1,1), linetype="dashed")

#为数据增加新的一列Label，将上调和下调差异表达前十的基因绘制在火山图中
#新加一列Label
deg.data$Label = ""

#对差异表达基因的pvalue值进行从小到大排序
deg.data <- deg.data[order(deg.data$pvalue),]

#高表达的基因中，选择pvalue最小的10个
up.genes <- head(deg.data$gene_symbol[which(deg.data$Group == "up-regulated")], 10)

#低表达的基因中，选择pvalue最小的10个
down.genes <- head(deg.data$gene_symbol[which(deg.data$Group == "down-regulated")], 10)

#将up.genes和down.genes合并
deg.top10.genes <- c(as.character(up.genes), as.character(down.genes))

#将top10.gens加入到Label中
deg.data$Label[match(deg.top10.genes, deg.data$gene_symbol)] <- deg.top10.genes

#参数说明：https://www.jianshu.com/p/674f90e020fa
#改变火山图点的颜色和坐标轴标注，使图片更美观
#绘制group1的最终火山图
#对输出的图保存至相应目录
pdf("火山图.pdf")
ggscatter(deg.data, x = "log2FoldChange", y = "logP", 
          color = "Group",
          palette = c("#2f5688", "#BBBBBB", "#CC0000"),
          size = 1,
          label =deg.data$Label,
          font.label = 8,
          repel =T,
          xlim = c(-15, 15),  # x坐标轴的范围
          xlab = "log2FoldChange",
          ylab = "-log10(P-value)",) + theme_base() +
  geom_hline(yintercept = 1.30, linetype="dashed") +
  geom_vline(xintercept = c(-1,1), linetype="dashed")
dev.off()

3. 对RNA_seq_2018_08_16进行GO和KEGG分析

# GO_KEGG -----------------------------------------------------------------

#参考文章：https://www.jianshu.com/p/435d863e0238，

#清空环境变量
rm(list=ls())

#安装包
#BiocManager::install("clusterProfiler")
#BiocManager::install("stringr")
#加载包
library(clusterProfiler)
library(stringr)
library(DOSE)
library(ggplot2)

#clusterProfiler 包里的一些默认作图方法，例如
#barplot(kegg)  #富集柱形图
#dotplot(kegg)  #富集气泡图
#cnetplot(kegg) #网络图展示富集功能和基因的包含关系
#emapplot(kegg) #网络图展示各富集功能之间共有基因关系
#heatplot(kegg) #热图展示富集功能和基因的包含关系
#Barplot画图参数详解:http://blog.sciencenet.cn/blog-1468811-939797.html

#参考物种的基因注释数据库：人类org.Hs.eg.db，果蝇org.Dm.eg.db，拟南芥org.At.tair.db，小鼠org.Mm.eg.db。
#下载参考人的基因注释库
#BiocManager::install("org.Hs.eg.db")
#加载小鼠的基因注释库
library(org.Hs.eg.db)

#准备输入数据
##待分析的数据就是一串基因名称，可以是ensembl_id、entrze_id或者symbol_id等类型
##读入差异基因的列表（此处根据dongfeng的实验结果，选取significant_gene[pvalue<0.05,abs(log2FC)>=1]进行作图分析
##对于上调、下调基因，需要手动分割成两个单独文件

#设置工作目录
setwd("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/")

#对样品进行读入
sig.gene_up <- read.csv("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/significant_pvalue_different_genes_genecount_up.csv")
sig.gene_dn <- read.csv("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/significant_pvalue_different_genes_genecount_down.csv")

##提取差异基因的列表
gene_up <- sig.gene_up$X
gene_dn <- sig.gene_dn$X

##调整数据格式为字符
gene_up <- as.character(gene_up)
gene_dn <- as.character(gene_dn)

##对基因由SYMBOL转换为ENTREZID格式
##select(x, keys, columns, keytype, ...)：基于keys, columns和keytype以data.frame数据类型返回数据，可以是一对多的关系
##mapIds(x, keys, column, keytype, ..., multiVals)： 类似于select，只不过就返回一个列。
gene_up.df <- select(org.Hs.eg.db,           #人基因注释库
                    keys = gene_up,          #样品组信息
                    columns = "ENTREZID",    #指定基因名称类型
                    keytype ="SYMBOL")       #输入的基因名称类型

gene_dn.df <- select(org.Hs.eg.db,           #小鼠基因注释库
                    keys=gene_dn,            #样品组信息
                    columns = "ENTREZID",    #指定基因名称类型
                    keytype="SYMBOL")        #输入的基因名称类型

#GO富集分析
GO_BP_up <- enrichGO(gene= gene_up.df$ENTREZID,          #基因列表文件中的基因名称
                     OrgDb = org.Hs.eg.db,               #指定物种的基因数据库，此为人
                     keyType = "ENTREZID",               #指定给定的基因名称类型，此为symbol_id
                     ont = "BP",                         #可选BP、MF、CC，也可以指定 ALL 同时计算 3 者
                     pAdjustMethod = "BH",               #指定p值校正方法
                     pvalueCutoff = 0.05,                #指定p值阈值，不显著的值将不显示在结果中
                     qvalueCutoff = 0.2,                 #指定q值阈值，不显著的值将不显示在结果中
                     readable = TRUE)                    #whether mapping gene ID to gene Name
# minGSSize                         minimal size of genes annotated by Ontology term for testing.
# maxGSSize                         maximal size of genes annotated for testing
# pool                              If ont=’ALL’, whether pool 3 GO sub-ontologies

GO_BP_dn <- enrichGO(gene= gene_dn.df$ENTREZID,          #基因列表文件中的基因名称
                     OrgDb = org.Hs.eg.db,               #指定物种的基因数据库，此为小鼠
                     keyType = "ENTREZID",               #指定给定的基因名称类型，此为symbol_id
                     ont = "BP",                         #可选BP、MF、CC，也可以指定 ALL 同时计算 3 者
                     pAdjustMethod = "BH",               #指定p值校正方法
                     pvalueCutoff = 0.05,                #指定p值阈值，不显著的值将不显示在结果中
                     qvalueCutoff = 0.2,                 #指定q值阈值，不显著的值将不显示在结果中
                     readable = TRUE)                    #whether mapping gene ID to gene Name 
# minGSSize                         minimal size of genes annotated by Ontology term for testing.
# maxGSSize                         maximal size of genes annotated for testing
# pool                              If ont=’ALL’, whether pool 3 GO sub-ontologies


#对GO富集分析的结果进行输出保存
#设置工作目录
setwd("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/GO_KEGG/")

write.csv(as.data.frame(GO_BP_up), "GO_BP_up.csv")
write.csv(as.data.frame(GO_BP_dn), "GO_BP_dn.csv")

#对GO富集分析进行绘图并输出保存
tiff("GO_BP_up.tiff")
barplot(GO_BP_up, showCategory = 16,title="The GO_BP enrichment analysis of all DEGs ") + 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(GO_BP_up) str_wrap(GO_BP_up, width = 30))
dev.off()

tiff("GO_BP_dn.tiff")
barplot(GO_BP_dn, showCategory = 16,title="The GO_BP enrichment analysis of all DEGs ") + 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(GO_BP_dn) str_wrap(GO_BP_dn, width = 30))
dev.off()

#KEGG富集分析
#clusterProfiler的KEGG富集分析方法特殊，它无需加载本地注释库，
#自动使用KEGG的在线数据库进行注释，因此给定的基因名称只能识别entrze id。
#每次打开R计算时，它会自动连接kegg官网获得最近的物种注释信息，因此数据库一定都是最新的

KEGG_up <- enrichKEGG( gene = gene_up.df$ENTREZID,            #基因列表文件中的基因名称
                       keyType = 'kegg',                      #kegg 富集
                       organism = 'hsa',                      #物种，mmu 代表小鼠，hsa代表人类，oas代表绵羊
                       pAdjustMethod = 'BH',                  #指定 p 值校正方法
                       pvalueCutoff = 0.05,                   #指定 p 值阈值，不显著的值将不显示在结果中
                       qvalueCutoff = 0.2,                    #指定 q 值阈值，不显著的值将不显示在结果中
                       use_internal_data= FALSE )

KEGG_dn <- enrichKEGG( gene = gene_dn.df$ENTREZID,            #基因列表文件中的基因名称
                       keyType = 'kegg',                      #kegg 富集
                       organism = 'hsa',                      #物种，mmu 代表小鼠，hsa代表人类，oas代表绵羊
                       pAdjustMethod = 'BH',                  #指定 p 值校正方法
                       pvalueCutoff = 0.05,                   #指定 p 值阈值，不显著的值将不显示在结果中
                       qvalueCutoff = 0.2,                    #指定 q 值阈值，不显著的值将不显示在结果中
                       use_internal_data= FALSE )

#KEGG分析结果各列内容：
#ID和Description，富集到的KEGG id和描述；
#GeneRatio和BgRatio，分别为富集到该KEGG条目中的基因数目/给定基因的总数目，以及该条目中背景基因总数目/该物种所有已知的KEGG功能基因数目；
#pvalue、p.adjust和qvalue，p值、校正后p值和q值信息；
#geneID和Count，富集到该KEGG条目中的基因名称（分析中使用的entrze id，故这里也显示的entrze id）和数目。
#如期望显示其它类型的基因id，如通俗的symbol id等类型，由于该分析中只能输入entrze id，因此可以通过基因名称转换的方式对entrze id和symbol id作个匹配转换。

#输出结果
write.table(KEGG_up, 'KEGG_up.csv', sep = ',', quote = FALSE, row.names = FALSE)
write.table(KEGG_dn, 'KEGG_dn.csv', sep = ',', quote = FALSE, row.names = FALSE)

#将输出结果的"ENTREZID"转为"SYMBOL"
#读取数据

#对group_1_up的gene进行名称类型转换
KEGG_up_df <- read.csv("KEGG_up.csv") 
#查看文件的整体信息以确认需要转换的列
KEGG_up_df
#with({})函数中的花括号语句，只针对括号内的语句执行，无需担心名字的冲突
#with()函数局限性在于赋值仅在此函数的括号内生效
with(KEGG_up_df, {
   
  gene = NA    #对变量gene进行一个定义
  for (i in 1:nrow(KEGG_up_df))        #for（i in seq)语句做循环，#nrow()对数据取行数
  {
   
    aaa = unlist(str_split(KEGG_up_df$geneID[i],"/"))  #str_split函数对字符串进行拆分，unlist函数将拆分的数据合并
    gene = c(gene, aaa)      #将所有的拆分结果合并在一起
  }
  KEGG_gene_up <<- gene[-1]        #对最开始的gene赋值进行剔除
})
#查看最终提取结果
KEGG_gene_up
#对提取的"ENTREZID"利用select()函数进行转换
#使用birt()函数也可以对ID和symbol进行转换，但是对于重复ID不再进行重复输出
KEGG_gene_up <- select(org.Hs.eg.db,               #人基因注释库
                       keys=KEGG_gene_up,          #样品组信息
                       columns = "SYMBOL",         #指定基因名称类型
                       keytype="ENTREZID")         #输入的基因名称类型
#查看最终转换结果
KEGG_gene_up

#对group_1_dn的gene进行名称类型转换
KEGG_dn_df <- read.csv("KEGG_dn.csv") 
#查看文件的整体信息以确认需要转换的列
KEGG_dn_df
#with({})函数中的花括号语句，只针对括号内的语句执行，无需担心名字的冲突
#with()函数局限性在于赋值仅在此函数的括号内生效
with(KEGG_dn_df, {
   
  gene = NA    #对变量gene进行一个定义
  for (i in 1:nrow(KEGG_dn_df))        #for（i in seq)语句做循环，#nrow()对数据取行数
  {
   
    aaa = unlist(str_split(KEGG_dn_df$geneID[i],"/"))  #str_split函数对字符串进行拆分，unlist函数将拆分的数据合并
    gene = c(gene, aaa)      #将所有的拆分结果合并在一起
  }
  KEGG_gene_dn <<- gene[-1]        #对最开始的gene赋值进行剔除
})
#查看最终提取结果
KEGG_gene_dn
#对提取的"ENTREZID"利用select()函数进行转换
#使用birt()函数也可以对ID和symbol进行转换，但是对于重复ID不再进行重复输出
KEGG_gene_dn <- select(org.Hs.eg.db,               #人基因注释库
                       keys=KEGG_gene_dn,          #样品组信息
                       columns = "SYMBOL",         #指定基因名称类型
                       keytype="ENTREZID")         #输入的基因名称类型
#查看最终转换结果
KEGG_gene_dn

#对转换后的结果输出保存
write.table(KEGG_gene_up, 'KEGG_gene_up.csv', sep = ',', quote = FALSE, row.names = FALSE)
write.table(KEGG_gene_dn, 'KEGG_gene_dn.csv', sep = ',', quote = FALSE, row.names = FALSE)

#对KEGG富集分析的结果进行输出保存
tiff("KEGG_up.tiff")
barplot(KEGG_up, showCategory = 16,title="The KEGG enrichment analysis of all DEGs ") + 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(KEGG_up) str_wrap(KEGG_up, width = 30))
dev.off()

tiff("KEGG_dn.tiff")
barplot(KEGG_dn, showCategory = 16,title="The KEGG enrichment analysis of all DEGs ") + 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(KEGG_dn) str_wrap(KEGG_dn, width = 30))
dev.off()

4.对RNA_seq_2018_08_16GO中分析上调的基因对应的promoter区进行motif分析

1. 根据gene_list从hg19_genes.gtf文件中提取gene_TSS信息

也可以使用biomaRt包进行注释，中间出错尚未解决，有机会再尝试。

# 提取待分析的gene_list
# 读入待处理文件
my_file <- read.csv("G:/dongfeng/RNA-seq/2018_08_16/Rtreatment/significant_different_genes/significant_pvalue_different_genes_genecount_up.csv")
# 根据需要提取gene_symbol
gene_symbol <- my_file$X

# 处理基因组注释文件
# 读入特定基因组文件
hg19_genes_gtf <- read.csv("G:/hg19_genes.gtf", header = FALSE, sep = "\t")
# 对读入的数据某一列进行拆分，此处以空格为例
library(stringr)
# 设置空的"hg19_genes_gtf_V1"向量，其行数与待处理数据行数一致
hg19_genes_gtf_V1<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V2<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V3<-rep(NA,nrow(hg19_genes_gtf))
hg19_genes_gtf_V4<-rep(NA,nrow(hg19_genes_gtf))
# 利用for循环，对hg19_genes_gtf数据框中的第9列进行拆分提取，每一个拆分元素储存至一个变量中
for (i in 1:nrow(hg19_genes_gtf)){
   
  hg19_genes_gtf_V1[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[1]
  hg19_genes_gtf_V2[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[2]
  hg19_genes_gtf_V3[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[3]
  hg19_genes_gtf_V4[i] <- unlist(str_split(hg19_genes_gtf[i,9], pattern = " "))[4]
}
# 对原数据框中的特定序列重新赋值
hg19_genes_gtf <- cbind(hg19_genes_gtf, hg19_genes_gtf_V1, hg19_genes_gtf_V2, hg19_genes_gtf_V3, hg19_genes_gtf_V4)

# 原始gtf文件中字符串末尾有分号故需要删除
# 此处以结尾最后一个字符为例，采用nchar()函数对字符串进行计数后减一，来代表最后一个字符。

# 设置空的"hg19_genes_gtf_V1"向量，其行数与待处理数据行数一致
# 此处NA可以用0代替
hg19_genes_gtf_V5<-rep(NA

最低0.47元/天解锁文章

垚垚爸爱学习

关注

1
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
操作记录-2020-12-21: dongfeng_GAS_project

今天老菜鸟接了实验室一个大活，要帮学生把一个课题结掉，因此需要处理三批次数据。当然了，老菜鸟这些天正在忙一些其他事情，不知道这些事情有了结果后，还会不会再继续这些工作了，暂且记录一下自己目前复杂的心情吧。1. 建立相应目录# 新建对应的目录mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt scripts_log2. 检查数据完整性这次将提供的md5.txt中数据提前编辑了一下，修改了样品
复制链接

扫一扫

专栏目录