拟南芥转录组分析

虽然现在已经不再做生信但是作为对之前学习的总结，记录一下转录组分析的流程。

 

流程参考陈铭老师的《生物信息学》第三版。

（一）数据预处理

1.提取fastq文件：使用fastq-dump工具从SRA文件中提取fastq文件

##使用fastq-dump将SRA文件转为fastq格式
fastq-dump -h
fastq-dump SRR3418005.sra &
fastq-dump SRR3418006.sra &
fastq-dump SRR3418019.sra &
fastq-dump SRR3418020.sra &

2.质量评估 使用FastQC检测原始测序数据质量

fastqc -o ../fastqc_results -f fastq SRR3418005.fastq 
SRR3418006.fastq SRR3418019.fastq SRR3418020.fastq &

3.质量控制 使用fast_trimmer截去reads前12位碱基，使用fastq_quality_filter筛选高质量reads

fastq_quality_filter -Q 33 -q 20 -p 80 -i SRR3418005_trimmed.fastq -o SRR3418005_filtered.fastq &
fastq_quality_filter -Q 33 -q 20 -p 80 -i SRR3418006_trimmed.fastq -o SRR3418006_filtered.fastq &
fastq_quality_filter -Q 33 -q 20 -p 80 -i SRR3418019_trimmed.fastq -o SRR3418019_filtered.fastq &
fastq_quality_filter -Q 33 -q 20 -p 80 -i SRR3418020_trimmed.fastq -o SRR3418020_filtered.fastq &

fastx_trimmer -Q 33 -f 12 -i SRR3418005.fastq -o fastx_results/ SRR3418005_trimmed.fastq &
fastx_trimmer -Q 33 -f 12 -i SRR3418006.fastq -o fastx_results/ SRR3418006_trimmed.fastq &
fastx_trimmer -Q 33 -f 12 -i SRR3418019.fastq -o fastx_results/ SRR3418019_trimmed.fastq &
fastx_trimmer -Q 33 -f 12 -i SRR3418020.fastq -o fastx_results/ SRR3418020_trimmed.fastq &

（二）比对和拼接

Hisat2 比对，首先用hisat2-build 构建index，然后比对，输出sam格式，再用samtools将sam转为bam格式，并排序构建索引。

1.建立索引并进行reads比对

hisat2-build ../genome/fasta/tair10.fasta tair10 #建立基因组索引库
hisat2 -p 2 --dta -x ../data/index/tair10 -U ../data/fastx_results/SRR3418005_filtered.fastq -S SRR3418005.sam &
hisat2 -p 2 --dta -x ../data/index/tair10 -U ../data/fastx_results/SRR3418006_filtered.fastq -S SRR3418006.sam &
hisat2 -p 2 --dta -x ../data/index/tair10 -U ../data/fastx_results/SRR3418019_filtered.fastq -S SRR3418019.sam &
hisat2 -p 2 --dta -x ../data/index/tair10 -U ../data/fastx_results/SRR3418020_filtered.fastq -S SRR3418020.sam

2.SAM文件处理

可见上一步生成的比对后的文件格式为.sam，下面就对sam文件进行操作。使用samtools对SAM文件排序并转换为BAM文件。

samtools sort -@ 2 -m 200M -o SRR3418005.bam SRR3418005.sam &
samtools sort -@ 2 -m 200M -o SRR3418006.bam SRR3418006.sam &
samtools sort -@ 2 -m 200M -o SRR3418019.bam SRR3418019.sam &
samtools sort -@ 2 -m 200M -o SRR3418020.bam SRR3418020.sam &

（三）计数

使用HTSeq-count从比对结果中提取所有基因匹配的reads数目

一般来说，RNA-seq不太会涉及到novel transcript的discovery，所以在hisat2比对完以后，都是衔接的reads counts提取，而不是stringtie（因为还没有成熟的流程配合）

htseq-count -q -f bam -s no -i gene_name ../../alignment/SRR3418005.bam ../../data/genome/gff/tair10.gtf > SRR3418005.count &
htseq-count -q -f bam -s no -i gene_name ../../alignment/SRR3418006.bam ../../data/genome/gff/tair10.gtf > SRR3418006.count &
htseq-count -q -f bam -s no -i gene_name ../../alignment/SRR3418019.bam ../../data/genome/gff/tair10.gtf > SRR3418019.count &
htseq-count -q -f bam -s no -i gene_name ../../alignment/SRR3418020.bam ../../data/genome/gff/tair10.gtf > SRR3418020.count &

四个文件合并成一个文本文件

join SRR3418005.count SRR3418006.count | join - SRR3418019.count | join - SRR3418020.count > count.txt #count文件整合
sed -i 's/ /\t/g' count.txt

（四）差异表达分析

这一步在Rstudio中实现

#差异表达分析


setwd("F:/rnaseq/data") #设置工作目录
library(DESeq2) #载入DESeq2包

#reads计数数据表操作
countTable <- read.table("count.txt", sep = "\t", header = FALSE) #读入reads计数矩阵
tail(countTable) #查看数据表最后的“小尾巴”
countTable <- countTable[- c(33611:33615), ] #去除描述行
rownames(countTable) <- countTable$V1 #将基因ID设置为行名
countTable <- countTable[, - 1] #删除基因ID列 
colnames(countTable) <- c("SRR3418005", "SRR3418006", "SRR3418019", "SRR3418020") #更改数据表列名
countTable <- countTable[- which(rowSums(countTable) < 4), ] #过滤count总数小于4的基因
nrow(countTable) #查看数据表行数（基因个数）

#设置样本处理信息（实验 vs. 对照）
colData <- data.frame(row.names = colnames(countTable), condition = c("ABA", "mock", "ABA", "mock"))

#DESeq2操作
dds <- DESeqDataSetFromMatrix(countData = countTable, colData = colData, design = ~ condition) #生成DESeqDataSet数据集WHICH
dds #查看数据集
dds$condition #查看样本处理信息
dds$condition <- relevel(dds$condition, "mock") #更改mock水平（使DESeq计算FoldChange时mock组作为分母）#Mock为对照组，默认A<M 所以要换mock为分母
dds$condition #查看更改水平后的样本处理信息

dds <- DESeq(dds) #差异表达计算
res <- results(dds) #生成差异表达结果表
summary(res) #查看总结信息（表达上调，下调等）
resOrdered <- res[order(res$padj), ] #按照校准后p值排序
write.csv(resOrdered, "DESeq2_results_all.csv") #将差异表达分析结果输出到csv文件

deg <- subset(resOrdered, padj <= 0.01 & abs(log2FoldChange) >= 2) #筛选显著差异表达基因（padj小于0.01且FoldChange绝对值大于4）（绝对值可以保证p>2和p<-2 两个方面）

summary(deg) #查看筛选后的总结信息
write.csv(deg, "DESeq2_results_significant.csv") #将差异表达显著的结果输出到csv文件

作图略