操作记录-2020-10-30：daizhongye_RNA_seq

最新推荐文章于 2021-05-29 10:59:54 发布

垚垚爸爱学习

最新推荐文章于 2021-05-29 10:59:54 发布

阅读量4.5k

点赞数 1

分类专栏：操作记录

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

版权

这篇博客详细记录了RNA-seq数据的分析流程，包括检查数据完整性，使用FastQC、Hisat2、SAMtools、Stringtie进行预处理，接着在Rstudio中使用DESeq2、ggpubr、pheatmap和clusterProfiler进行差异基因分析、绘制火山图、heatmap图及GO_KEGG富集分析。

摘要由CSDN通过智能技术生成

隔了一个月又来操作一次，老菜鸟的脑子明显不够用啊。

dir=/f/xudonglab/zexing/projects/daizhongye/RNA_seq/2020_10_29
mkdir -p ${dir}/aligned ${dir}/aligned/ballgown ${dir}/aligned/bam ${dir}/aligned/bam.index ${dir}/aligned/bam.sort ${dir}/aligned/sam ${dir}/fastqc_report ${dir}/GSEA ${dir}/MD5_txt ${dir}/raw ${dir}/scripts_log

1. 检查上传数据的完整性

操作记录如下：

#显示各数据的完整性代码
(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/MD5_txt$ cat md5.txt
efcc57e9ff0bd5c3f4a04dcd01513903  raw/E14_Scr_SL_1.fq.gz
657fe81df85e0f3c246c3bcfbc07a4c3  raw/E14_Scr_SL_2.fq.gz
56808aabfa366957e63df968e586d7bc  raw/E14_shT1_1_1.fq.gz
6497464b97496d6da76d22177bddd672  raw/E14_shT1_1_2.fq.gz
e292972349e7331dc28200421b6a927c  raw/E14_shT1_2_1.fq.gz
366826927ab1ea176d82352daf08f6ba  raw/E14_shT1_2_2.fq.gz
0f132bcf9b81fecb450db1ba33b2f577  raw/E14_shT2_1_1.fq.gz
b16d41d1eabcb03feafce2fd197090b8  raw/E14_shT2_1_2.fq.gz
f201561de55191ebfcf8987f1ae45571  raw/E14_shT2_2_1.fq.gz
99693aa2350c21c1f57b9845740c604e  raw/E14_shT2_2_2.fq.gz
#将各数据的完整代码统一写入待检查的文件中
(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/MD5_txt$ echo "efcc57e9ff0bd5c3f4a04dcd01513903  raw/E14_Scr_SL_1.fq.gz
657fe81df85e0f3c246c3bcfbc07a4c3  raw/E14_Scr_SL_2.fq.gz
56808aabfa366957e63df968e586d7bc  raw/E14_shT1_1_1.fq.gz
6497464b97496d6da76d22177bddd672  raw/E14_shT1_1_2.fq.gz
e292972349e7331dc28200421b6a927c  raw/E14_shT1_2_1.fq.gz
366826927ab1ea176d82352daf08f6ba  raw/E14_shT1_2_2.fq.gz
0f132bcf9b81fecb450db1ba33b2f577  raw/E14_shT2_1_1.fq.gz
b16d41d1eabcb03feafce2fd197090b8  raw/E14_shT2_1_2.fq.gz
f201561de55191ebfcf8987f1ae45571  raw/E14_shT2_2_1.fq.gz
99693aa2350c21c1f57b9845740c604e  raw/E14_shT2_2_2.fq.gz" >check_md5sum.txt
#执行md5sum 命令对其完整性进行检测
(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29$ md5sum -c check_md5sum.txt
raw/E14_Scr_SL_1.fq.gz: OK
raw/E14_Scr_SL_2.fq.gz: OK
raw/E14_shT1_1_1.fq.gz: OK
raw/E14_shT1_1_2.fq.gz: OK
raw/E14_shT1_2_1.fq.gz: OK
raw/E14_shT1_2_2.fq.gz: OK
raw/E14_shT2_1_1.fq.gz: OK
raw/E14_shT2_1_2.fq.gz: OK
raw/E14_shT2_2_1.fq.gz: OK
raw/E14_shT2_2_2.fq.gz: OK

2. 使用FastQC软件对数据进行质控

vim新建fastqc_script脚本如下：

#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for analysis of RNA-seq data by FastQC.
#History:
# 2020/10/29         zexing              First release
#fastqc命令为质控命令
#Usage: fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
#简写代码：fastqc -t 16 -o <out-dir> seqfile1
#调用程序fastqc，参数-t设置线程数为8，参数-o设置结果输出的目录，参数-c可以加入污染物选项（接头信息），最后为读入文件。
dir=/f/xudonglab/zexing/projects/daizhongye/RNA_seq/2020_10_29
fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/raw/*.fq.gz

后台运行fastqc_scrip如下脚本：

(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/scripts_log$ nohup bash fastqc_script > fastqc_script_log &
[1] 176237

3. 使用Hisat2软件对测序结果进行比对

vim新建hisat2_script脚本如下：

#! /bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for aligning of RNA-seq data.
#History:
# 2020/10/30         zexing              First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#hisat2命令为序列比对命令
#Usage: hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]
#简写代码：hisat2 -t -p 8 -x <ht2-idx> -1 fq.gz -2 fq.gz -S <sam>
#调用程序hisat2，参数-t显示时间，参数-p设定线程数，参数-x指定参考基因组索引文件的前缀（目录及文件前缀）
#参数-1指定双向测序的第一条.fq.gz测序数据，如果多组数据，使用逗号将文件分隔
#参数-2指定双向测序的第二条.fq.gz测序数据，如果多组数据，使用逗号将文件分隔
#参数-S指定比对结果的输出目录及名称
dir=/f/xudonglab/zexing/projects/daizhongye/RNA_seq/2020_10_29
for i in E14_Scr_SL E14_shT1_1 E14_shT1_2 E14_shT2_1 E14_shT2_2
do
hisat2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_mm10/hisat2_index/hisat2_index_mm10 \
-1 ${dir}/raw/${i}_1.fq.gz \
-2 ${dir}/raw/${i}_2.fq.gz \
-S ${dir}/aligned/sam/${i}.sam
done

后台运行hisat2_script脚本如下：

(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/scripts_log$ nohup bash hisat2_script > hisat2_script_log &
[1] 180160

4. 使用SAMtools软件对文件进行格式转换、排序

vim新建samtools_script如下：

#!/bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#这里写一个脚本，将samtools的view、sort、index和flagstat四个命令串联在一起，对RNA-seq数据批量处理。
#History:
# 2020/10/30         李泽兴      First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#samtools view为格式转换命令
#Usage: samtools view [options] <in.bam> | <in.sam> | <in.cram> [region...]
#-@设置线程数为8， -S输入sam文件，-1b使用快速压缩比生成bam文件，-o设定标准输出文件名。
#简写代码：samtools view -@ 8 -S {$i}.sam -1b -o {$i}.bam
#samtools sort为排序命令
#Usage: samtools sort [options..] [in.bam]
#-@设置线程数为8，-l设置压缩比为5，-o设定标准输出文件名，最后一个为输入.bam文件，该命令默认为染色体位置排序，不影响后期操作。
#简写代码：samtools sort -@ 8 -l 5 -o {$i}.bam.sort {$i}.bam
#index为建立索引命令，目前还不知道这样建立索引有何用处，暂且写下来以备后续使用。
#Usage: samtools index [-bc] [-m INT] <in.bam> <out.index>
#-@设置线程数为8，<in.bam>输入bam文件，<out.index>输出index文件。
#简写代码：samtools index -@ 8 {$i}.bam.sort {$i}.bam.index
#flagstat为查看比对情况命令，目前暂且对这部分不是很熟悉，暂且写下来以备后续使用。
#Usage: samtools flagstat [options] <in.bam>
#-@设置线程数为8，对于输出内容重定向写入>指定文件
#简写代码：samtools flagstat -@ 8 {$i}.bam.sort > {$i}.bam.flagstat
dir=/f/xudonglab/zexing/projects/daizhongye/RNA_seq/2020_10_29/aligned
for i in E14_Scr_SL E14_shT1_1 E14_shT1_2 E14_shT2_1 E14_shT2_2
do
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 ${dir}/bam.index/${i}.bam.index
samtools flagstat -@ 16 ${dir}/bam.sort/${i}.bam.sort > ${dir}/${i}.bam.flagstat
done

后台运行samtools_script脚本如下：

#后台运行命令如下：
 zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/scripts_log$ nohup bash samtools_script > samtools_script_log &
[3] 6353

5. 使用Stringtie软件进行拼接和定量

vim新建stringtie_script如下：

#!/bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
#program
#     This program is to perform HTSeq-count assay for RNA-seq data.
#History
#     2020/10/30       zexing            First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#创建针对不同样品的ballgown/"$i"目录，以存放输出结果
#stringtie为转录本拼接和定量的软件
#利用stringtie分别输出gene count数目和FPKM值，gene count数用户基因差异分析，FPKM值为基因表达值，用于heatmap等绘制
#Usage: stringtie <aligned_reads.bam> [options]
#aligned_reads.bam 是输入文件，该输入文件要求必须按基因组位置排序
#参数-o [<path/>]<out.gtf>，设置StringTie组装转录本的输出GTF文件的路径和文件名。
#参数-p <int>指定组装转录本的线程数（CPU）
#参数-G <ref_ann.gff>使用参考注释基因文件指导组装过程，格式GTF/GFF3。输出文件中既包含已知表达的转录本，也包含新的转录本。
#参数-B 应用该选项，则会输出Ballgown输入表文件（* .ctab），其中包含用-G选项给出的参考转录本的覆盖率数据。
#参数-e 限制reads比对的处理，仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项，则会跳过处理与参考转录本不匹配的组装转录本，这将大大的提升了处理速度。
#参数-A <gene_abund.tab> 输出基因丰度的文件（制表符分隔格式）
#简写代码：stringtie {$i}.bam.sort -o path/name.gtf -p 16 -G gene.gtf -eB -A path/name.tab
#此次主要输出两个文件：gtf文件用于转换gene count文件；tab文件包含FPKM数值。
dir=/f/xudonglab/zexing/projects/daizhongye/RNA_seq/2020_10_29/aligned
for i in E14_Scr_SL E14_shT1_1 E14_shT1_2 E14_shT2_1 E14_shT2_2
do
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_mm10/mm10_genes.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done

后台运行stringtie_script脚本如下

(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/scripts_log$ nohup bash stringtie_script > stringtie_script_log &
[1] 36868

6. 使用prepDE.py脚本提取read counts数值

注意：prepDE.py为python2的脚本，应该先将python设置为低版本后再运行脚本。
操作记录如下：

#从之前的文件夹中将prepDE.py脚本拷贝至当前文件夹中的ballgown子文件夹中
(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/aligned/ballgown$  cp /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/ballgown_1/prepDE.py ./

#退出当前conda环境
(base) zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/aligned/ballgown$ conda deactivate

#检查服务器中的python版本信息
zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/aligned/ballgown$ python --version
Python 2.7.12

#使用python命令直接运行脚本
zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/aligned/ballgown$ python prepDE.py

#查看运行结果
zexing@DNA:~/projects/daizhongye/RNA_seq/2020_10_29/aligned/ballgown$ ll
total 1.8M
drwxrwxr-x 7 zexing zexing 4.0K 10月 30 13:02 .
drwxrwxr-x 8 zexing zexing 4.0K 10月 29 16:26 ..
drwxrwxr-x 3 zexing zexing 4.0K 10月 30 12:46 E14_Scr_SL
drwxrwxr-x 2 zexing zexing 4.0K 10月 30 12:48 E14_shT1_1
drwxrwxr-x 2 zexing zexing 4.0K 10月 30 12:51 E14_shT1_2
drwxrwxr-x 2 zexing zexing 4.0K 10月 30 12:55 E14_shT2_1
drwxrwxr-x 2 zexing zexing 4.0K 10月 30 12:59 E14_shT2_2
-rw-rw-r-- 1 zexing zexing 794K 10月 30 13:02 gene_count_matrix.csv
-rw-rw-r-- 1 zexing zexing  12K 10月 30 12:47 prepDE.py
-rw-rw-r-- 1 zexing zexing 987K 10月 30 13:02 transcript_count_matrix.csv

#其中"gene_count_matrix.csv"即是DESeq2的输入文件。

以下分析在本地机的Rstudio中完成

7. 利用DESeq2包进行差异基因分析

代码如下：

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

#清空环境变量
rm(list=ls())
#设置工作目录
setwd("G:/daizhongye/RNA-seq/2020_10_29/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)
#将各组数据分开
gene_count_group_1 <- gene_count[, 1:3]
gene_count_group_2 <- gene_count[, c(1, 4, 5)]
#将该文件保存至对应目录
write.csv(gene_count_group_1, file = "G:/daizhongye/RNA-seq/2020_10_29/gene_count/gene_count_group_1.csv", row.names = TRUE)
write.csv(gene_count_group_2, file = "G:/daizhongye/RNA-seq/2020_10_29/gene_count/gene_count_group_2.csv", row.names = TRUE)

#加载DESeq2包
library(DESeq2)
#DESeq2需要三种矩阵，分别为countData表达矩阵,colData样品信息矩阵及design差异表达矩阵
#countData为表达矩阵即gene_count
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#设置condition样品组别、重复数
condition_group_1 <- factor(c(rep("scr", 1), rep("shT1", 2)), levels = c("scr","shT1"))
condition_group_2 <- factor(c(rep("scr", 1), rep("shT2", 2)), levels = c("scr","shT2"))
#显示condition设置信息
condition_group_1
condition_group_2
#设置group组对应的样品信息矩阵colData
colData_group_1 <- data.frame(row.names = colnames(gene_count_group_1), condition_group_1)
colData_group_2 <- data.frame(row.names = colnames(gene_count_group_2), condition_group_2)
#显示colData设置信息
colData_group_1
colData_group_2
#在R里面用于构建公式对象，~左边为因变量，右边为自变量。
#标准流程：dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) 
#countData为表达矩阵即countdata
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#对dds_group进行标准流程构建
dds_group_1 <- DESeqDataSetFromMatrix(gene_count_group_1, colData_group_1, design = ~condition_group_1)
dds_group_2 <- DESeqDataSetFromMatrix(gene_count_group_2, colData_group_2, design = ~condition_group_2)
#对原始dds_group进行normalize
dds_group_1 <- DESeq(dds_group_1)
dds_group_2 <- DESeq(dds_group_2)
#显示dds信息
dds_group_1
dds_group_2

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

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

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

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

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

#将差异分析的所有结果进行输出保存
write.csv(res_group_1, file="G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/all_different_genes/all_different_genes_group_1_genecount.csv")
write.csv(res_group_2, file="G:/daizhongye/RNA-seq/2020_10_29/Rtreatment/all_different_genes/all_different_genes_group_2_genecount.csv")

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