RNA-seq分析（STAR-RSEM-DESeq2）
2021-12-13 本文已影响0人
清零2000
一、准备RNA-seq数据
- 需要参考基因组文件，注释文件，样本的fastq.gz文件
- 检查fastq.gz是否完整
md5sum sample.fastq.gz > md5.txt
md5sum -c md5.txt
#结果为OK既可，否则再次准备RNA-seq文件
二、fastqc生成质量报告（multiqc 提高查看效率）
- fastqc为每个样本生成质量报告，可视化的报告测序数据的质量，为下一步质量控制做准备。
- multiqc可以将多个fastqc得到的质量报告合成一份质量报告。如果有很多fastqc文件时，一份一份翻看效率很低，此时用mutilqc合并成一份报告可以提高效率
1. 安装fastqc和multiqc
conda install fastqc
conda install multiqc
2. 使用fastqc生成质量报告
##fastqc的几种使用方式
#多线程同时处理
ls */*fastq.gz | xargs -I [] nohup fastqc [] &
#批量处理
fastqc *fastq.gz
#单个处理
fastqc Sample.fastq.gz
3.使用multiqc命令合并指定目录下的html文件，并生成一个html文件(multiqc可以自动检测目录下的html文件，所以只需输入目录）
multiqc <指定目录>
三、trimmomatic质量控制
- 顾名思义，根据fastqc生成的质量报告对测序数据进行质量控制。比如，去除接头和低质量碱基，以减小噪点对比对和定量的影响
1.安装trimmomatic
conda install trimmomatic
2.trimmomatic参数
PE/SE 设定对Paired-End或Single-End的reads进行处理，其输入和输出参数稍有不一样。
-threads 设置多线程运行数
-phred33 设置碱基的质量格式，可选pred64，不设置这个参数，软件会自动判断输入文件是哪种格式
-trimlog 指定过滤日志文件名，日志中包含以下四列内容:read ID、过滤之后剩余序列长度、过滤之后的序列起始碱基位置（序列开头处被切掉了多少个碱基）、过滤之后的序列末端碱基位置、序列末端处被剪切掉的碱基数。
#由于生成的trimlog文件中包含了每一条 reads 的处理记录，因此文件体积巨大（GB级别），如果后面不会用到 trim日志，建议不要使用这个参数
-basein 通常 PE 测序的两个文件，指定其中 R1 文件名，软件会推测出 R2 的文件名，但是这个功能实测并不好用，建议不用-basein参数，直接指定两个文件名（R1 和 R2）作为输入
-baseout 输出文件有四个，使用 -baseout 参数指定输出文件的 basename，软件会自动为四个输出文件命名，过滤之后双端序列都保留的就是 paired，反之如果其中一端序列过滤之后被丢弃了另一端序列保留下来了就是 unpaired(即 成对的clean reads, 未成对的正向序列以及未成对的反向序列)
#一般情况下，若paired reads百分比占90%以上，可只对paired reads进行比对分析
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true
切除adapter序列。参数后面分别接adapter序列的fasta文件：第一步 seed 搜索时允许的最大错配碱基个数2：palindrome模式下匹配碱基数阈值30：simple模式下的匹配碱基数阈值10(7-15之间）：palindrome 模式允许切除的最短接头序列为 8bp（默认值）：palindrome 模式去除与 R1 完全反向互补的 R2（默认去除false），但在有些情况下，例如需要用到 paired reads 的 bowtie2 流程，就要将这个参数改为 true，否则会损失一部分 paired reads。
#按照规定顺序，ILLUMINACLIP 各个参数之间以冒号分开，PE测序需要注意最后一个参数。对于SE测序最后两个参数可以不设置
LEADING:3 切除首端碱基质量小于3的碱基
#Illumina平台有些低质量的碱基质量值被标记为 2 ，所以设置为3可以过滤掉这部分低质量碱基。
TRAILING:3 切除尾端碱基质量小于3的碱基
SLIDINGWINDOW:15:20
滑窗质量过滤，一般一个read的低质量序列都是集中在末端，也有很少部分在开头。从5'端开始进行滑动，当滑动位点周围一段序列(window)的平均碱基低于阈值，则从该处进行切除。Windows的size是15个碱基（一般设置在10-30之间），其平均碱基质量小于20，则切除
MINLEN:50 可被保留的最短reads长度，应根据原始序列的长度而定
HEADCROP:<length> 在reads的首端切除指定的长度
CROP:<length> 保留reads到指定的长度
TOPHRED33 将碱基质量转换为pred33格式
TOPHRED64 将碱基质量转换为pred64格式
3.使用trimmomatic质量控制
trimmomatic PE -threads 16 \ #PE指双端测序 -threads指使用几个线程
./01raw_data/sample_R1.fastq.gz ./01raw_data/sample_R2.fastq.gz \ #待质控样本的reads1和reads2文件
./02clean_data/sample_paired_clean_R1.fq.gz \#reads1中可以和read2配对的reads
./02clean_data/sample_unpair_clean_R1.fq.gz \#reads1中不可以和read2配对的reads
./02clean_data/sample_paired_clean_R2.fq.gz \#reads2中可以和read1配对的reads
./02clean_data/sample_unpair_clean_R2.fq.gz \#reads2中不可以和read1配对的reads
ILLUMINACLIP:~/miniconda3/share/trimmomatic-0.39-2/adapters/TruSeq2-PE.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
4.可能出现的接头不适配的问题
①判断需要哪种接头文件
fastqc报告会给出可能的接头来源，如下图
接头来源
- "Possible Source"是"Illumina Single End" 或 "Illumina Paired End" 时，根据自己是双端测序还是单端测序选择"TruSeq2-SE.fa" 或 "TruSeq2-PE.fa"
- "Possible Source"是"TruSeq Universal Adapter" 或 "TruSeq Adapter, index..." 时，根据自己是双端测序还是单端测序选择"TruSeq3-SE.fa" 或 "TruSeq3-PE.fa"
②接头文件选择
(base) [xxx@ln01 adapters]$ pwd
#接头文件路径
/home/xxx/miniconda3/share/trimmomatic-0.39-2/adapters
(base) [xxx@ln01 adapters]$ ls -lh
#接头文件一共有6类，根据fastqc给的接头报告选择相应的接头文件
-rw-rw-r-- 2 xxx xxx 239 Mar 28 2021 NexteraPE-PE.fa
-rw-rw-r-- 2 xxx xxx 538 Mar 28 2021 TruSeq2-PE.fa
-rw-rw-r-- 2 xxx xxx 142 Mar 28 2021 TruSeq2-SE.fa
-rw-rw-r-- 2 xxx xxx 259 Mar 28 2021 TruSeq3-PE-2.fa
-rw-rw-r-- 2 xxx xxx 93 Mar 28 2021 TruSeq3-PE.fa
-rw-rw-r-- 2 xxx xxx 119 Mar 28 2021 TruSeq3-SE.fa
- 把所需的接头文件的绝对路径放在 "ILLUMINACLIP:" 之后
四、STAR比对
- STAR可以将质控后的fastq.gz文件比对到参考基因组上，为后续的定量准备文件。STAR比对时，需要先用STAR构建他自己可以识别的索引文件，之后再比对
1.安装STAR
conda install STAR
2.构建索引
STAR --runThreadN 8 \#使用的线程，可根据电脑配置自行更改
--runMode genomeGenerate \#工作模式（建立索引）
--genomeDir 03star_index/ \#索引文件存储位置
--genomeSAindexNbases 12 \#报错，14太大，需要用12
--genomeFastaFiles 00ref/Brapa.genome.fasta \ #参考基因组
--sjdbGTFfile 00ref/Brapa.genome.gtf \ #注释文件
--sjdbOverhang 149 #一般为reads长度-1（150-1）
3.比对
STAR --runThreadN 16 --genomeDir 03star_index/ \ #索引文件所在目录
--readFilesCommand zcat \ #执行解压缩
--readFilesIn02clean_data/sample_paired_clean_R1.fastq.gz \
02clean_data/sample_paired_clean_R2.fastq.gz \ #输入文件
--outFileNamePrefix 04STAR_out/sample \#输出目录及文件前缀
--outSAMtype BAM SortedByCoordinate \#指定输出bam 文件并排序
--outBAMsortingThreadN 16 \ #设置排序中的线程数
--quantMode TranscriptomeSAM GeneCounts #将上面基因组比对的信息转化为转录本比对的信息,count
# --quantMode TranscriptomeSAM will output alignments translated into tran-script coordinates,为了使用RSEM 进行定量分析做准备
4.输出文件
- 一共有7个文件
-rw-r--r-- 1 xxx xxx 2.7G Dec 5 00:43 sampleAligned.sortedByCoord.out.bam
-rw-r--r-- 1 xxx xxx 2.4G Dec 5 00:34 sampleAligned.toTranscriptome.out.bam
-rw-r--r-- 1 xxx xxx 2.0K Dec 5 00:43 sampleLog.final.out
-rw-r--r-- 1 xxx xxx 1.4M Dec 5 00:43 sampleLog.out
-rw-r--r-- 1 xxx xxx 6.5K Dec 5 00:43 sampleLog.progress.out
-rw-r--r-- 1 xxx xxx 789K Dec 5 00:35 sampleReadsPerGene.out.tab
-rw-r--r-- 1 xxx xxx 6.1M Dec 5 00:35 sampleSJ.out.tab
五、RSEM定量
- RSEM定量时，同样需要先构建RSEM可以识别的索引文件，之后才可以进行定量
1.安装
conda install rsem
2.构建索引
rsem-prepare-reference \
--gtf 00ref/Brapa.genome.gtf \#参考基因组的注释文件
00ref/Brapa.genome.fasta \#参考基因组
05rsem_index/Brapa_genome #索引文件存储位置
3.定量
rsem-calculate-expression \
--paired-end \ #双端测序
--no-bam-output \ #不输出bam结果(bam是一种文件格式，sam -> bam）
--alignments -p 16 -q ./04STAR_out/sampleAligned.toTranscriptome.out.bam \#比对所用的线程和输入文件
./05rsem_index/Brapa_genome \#索引文件所在目录及其前缀
./06rsem_out/sample_rsem #输出目录及文件前缀
4.输出文件
- 一共有2个文件，一个目录。其中的两个文件sample_rsem.genes.results 和 sample_rsem.isoforms.results分别记录gene水平和转录本水平的定量结果。区别在于gene_id和transcript_id在第一列还是第二列
-rw-r--r-- 1 xxx xxx 2.2M Dec 5 11:39 sample_rsem.genes.results
-rw-r--r-- 1 xxx xxx 2.4M Dec 5 11:39 sample_rsem.isoforms.results
drwxr-xr-x 1 xxx xxx 4.0K Dec 5 11:39 sample_rsem.stat
- sample_rsem.genes.results
gene_id transcript_id(s) length effective_length expected_count TPM FPKM
Bra000001 Bra000001.1 1515.00 1319.96 2.00 0.08 0.10
Bra000002 Bra000002.1 1089.00 893.96 45.00 2.63 3.41
Bra000003 Bra000003.1 570.00 375.00 115.00 16.05 20.77
Bra000004 Bra000004.1 375.00 180.60 42.00 12.17 15.75
Bra000005 Bra000005.1 1434.00 1238.96 1522.50 64.29 83.25
Bra000006 Bra000006.1 312.00 118.83 613.00 269.90 349.45
Bra000007 Bra000007.1 711.00 515.97 59.00 5.98 7.75
Bra000008 Bra000008.1 882.00 686.96 12.00 0.91 1.18
Bra000009 Bra000009.1 657.00 461.97 1478.77 167.48 216.85
- sample_rsem.isoforms.results
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct
Bra000001.1 Bra000001 1515 1319.96 2.00 0.08 0.10 100.00
Bra000002.1 Bra000002 1089 893.96 45.00 2.63 3.41 100.00
Bra000003.1 Bra000003 570 375.00 115.00 16.05 20.77 100.00
Bra000004.1 Bra000004 375 180.60 42.00 12.17 15.75 100.00
Bra000005.1 Bra000005 1434 1238.96 1522.50 64.29 83.25 100.00
Bra000006.1 Bra000006 312 118.83 613.00 269.90 349.45 100.00
Bra000007.1 Bra000007 711 515.97 59.00 5.98 7.75 100.00
Bra000008.1 Bra000008 882 686.96 12.00 0.91 1.18 100.00
Bra000009.1 Bra000009 657 461.97 1478.77 167.48 216.85 100.00
六、构建表达矩阵
可以直接使用RSEM构建表达矩阵。表达矩阵是以sample_rsem.genes.results的expected_count列构建的
1.构建矩阵
rsem-generate-data-matrix *rsem.genes.results > output. matrix
2.输出结果
(base) [xxx@ln01 L7_24h-vs-L7_CK]$ head L7_24h-vs-L7_CK.matrix
"16group/L7_24h-vs-L7_CK/L7_24h_1/L7_24h_1rsem.genes.results" "16group/L7_24h-vs-L7_CK/L7_24h_2/L7_24h_2rsem.genes.results" "16group/L7_24h-vs-L7_CK/L7_24h_3/L7_24h_3rsem.genes.results" "16group/L7_24h-vs-L7_CK/L7_CK_1/L7_CK_1rsem.genes.results" "16group/L7_24h-vs-L7_CK/L7_CK_2/L7_CK_2rsem.genes.results" "16group/L7_24h-vs-L7_CK/L7_CK_3/L7_CK_3rsem.genes.results"
"Bra000001" 2.00 2.00 8.00 8.00 19.00 21.00
"Bra000002" 38.00 45.00 46.00 67.00 62.00 46.00
"Bra000003" 169.00 115.00 93.00 185.00 161.00 183.00
"Bra000004" 17.00 42.00 16.00 22.00 21.00 26.00
"Bra000005" 1425.79 1522.50 1186.24 1151.77 1059.85 1278.12
"Bra000006" 753.00 613.00 567.00 525.00 439.00 619.00
"Bra000007" 35.00 59.00 33.00 55.00 57.00 42.00
"Bra000008" 1.00 12.00 8.00 12.00 10.00 16.00
"Bra000009" 1651.77 1478.77 1346.60 666.35 831.43 959.23
表达矩阵的列名是以输入文件的路径所命名的，所以需要根据自己的需求进行修改，提供一个我修改列名的代码
sample="L7_24h-vs-L7_CK"
awk 'BEGIN{printf"geneid\t'${sample%%-*}'_1\t'${sample%%-*}'_2\t'${sample%%-*}'_3\t'${sample##*-}'_1\t'${sample##*-}'_2\t'${sample##*-}'_3\n"} /Bra/{print $0}' ./16group/${sample}/${sample}.matrix \
> ./16group/${sample}/${sample}_deseq2.txt
七、差异基因分析
1.在Rstudio中安装DESeq2
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
library(BiocManager)
install.packages("DESeq2")
2.分析
##批量处理
library(DESeq2)
library(stringr)
library(dbplyr)
file_name <- c("L7_24h-vs-L7_CK","L7_24h-vs-L7_Re24h",
"L7_24h-vs-Le_24h","L7_3h-vs-L7_24h",
"L7_3h-vs-L7_CK","L7_3h-vs-L7_Re24h",
"L7_3h-vs-Le_3h","L7_CK-vs-Le_CK",
"L7_Re24h-vs-L7_CK","L7_Re24h-vs-Le_Re24h",
"Le_24h-vs-Le_CK","Le_24h-vs-Le_Re24h",
"Le_3h-vs-Le_24h","Le_3h-vs-Le_CK",
"Le_3h-vs-Le_Re24h","Le_Re24h-vs-Le_CK")
for ( i in 1:16 ) {
##提取所比较的样本名字，在控制条件部分设置level。level高的为控制组，低的为参照组
#将 - 看为分隔符，将file_name分为三个部分，[[1]][1]为第一部分，[[1]][3]为第三部分
sample <- str_split(file_name[i],"-")
sample1 <- sample[[1]][1]
sample2 <- sample[[1]][3]
# read.table ,文件有header,第一列为row.name
input_data <- read.table(paste0(file_name[i],'_deseq2.txt'),
header = TRUE, row.names = 1)
# 取整,函数round
round_data <-round(input_data,digits = 0)
# 准备文件
# as.matrix 将输入文件转换为表达矩阵
count_data <- as.matrix(round_data)
# 控制条件：因子。level高的为控制组，低的为参照组。logFC=(控制组) / (参照组)
condition <- factor(c(rep(sample1, 3), rep(sample2, 3)),levels = c(sample2,sample1))
# input_data根据控制条件构建data.frame
coldata <- data.frame(row.names=colnames(count_data),
condition)
# build deseq input matrix 构建输入矩阵
#countData作为矩阵的input_data；colData Data.frame格式；
#控制条件design;
dds <- DESeqDataSetFromMatrix(
countData=count_data,colData=coldata, design=~condition)
#构建dds对象
dds <- DESeq(dds)
# 提取结果
#dds dataset格式转换为DESeq2 中result数据格式，矫正值默认0.1
res <- results(dds,alpha=0.05)
#查看res(DESeqresults格式),可以看到上下调基因
summary(res)
# 将经过矫正的表达量结果加进去；
res$up_down <- ifelse(res$log2FoldChange=="NA", "NA", ifelse(res$log2FoldChange<0, "down", "up"))
#resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
resdata <- merge(as.data.frame(res),as.data.frame(counts(dds,normalized=TRUE)),
by="row.names",sort=FALSE)
#将fpkm加入分析结果中
#resdata <- merge(as.data.frame(res), as.data.frame(read.table(paste0(file_name[i],'_fpkm.txt'),
# header = TRUE, row.names = #1)),by="row.names",sort=FALSE)
names(resdata)[1] <- "Gene"
# 输出结果,保存为csv文件
resdata$padj[is.na(resdata$padj)] <- 1
write.csv(resdata,file=paste0(file_name[i],'.csv'),sep = "\t",quote = F, row.names = F)
#筛选出 padj < 0.05, logFC >= 1 的基因
deg <- subset(resdata, resdata$padj <= 0.05 & abs(resdata$log2FoldChange) >= 1 )
#当与fpkm合并后，可能会出现fpkm的和为零但是显著差异的情况，这种基因需要去掉。
#deg <- subset(resdata, resdata$padj <= 0.05 & abs(resdata$log2FoldChange) >= 1 & apply(resdata[,9:14], 1, sum) > 0 )
write.csv(deg,file=paste0(file_name[i],'-deg-padj-0.05-logFC-1.csv'),sep = "\t",quote = F, row.names = F)
}
