RNA-seq数据分析

RNA-seq hisat2 —— stringtie —— DESeq2 pipeline

1、需求

1.1 软件

linux安装：
hisat2 conda install hisat2
stringtie conda install stringtie
samtools conda install samtools=1.15.1
gffread conda install gffread

R包：
DESeq2 BiocManager::install('DESeq2')

1.2 数据

(1) fastq测序数据

如图中例子，有9个样本，每个样本3个重复，双端测序

数据质控略

(2) 基因组及其gff注释文件，下载自Phytozome
gff转为gtf备用： gffread Wm82.gff3 -T -o Wm82.gtf

2、比对

2.1 构建索引

hisat2-build -p 8 path/to/genome.fa Wm82

2.2 准备config文件

ls path/to/fastq/*R1* > 1
ls path/to/fastq/*R2* > 2
ls path/to/fastq/*R2* | cut -d '/' -f 6 | cut -d '_' -f 1-3 > 0
paste 0 1 2 > align.config

2.3 比对

#!/bin/bash
hisat_index=/pub4/zhangchao/Methylome_Soybean/Reference/index/Wm82
cat align.config | while read id;
do
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
hisat2 -p 20 -x ${hisat_index} -1 $fq1 -2 $fq2 -S ${sample}.sam 2>${sample}_hisat.log
grep "NH:i:1" ${sample}.sam > ${sample}_unique.sam      #仅保留唯一比对
samtools view -H ${sample}.sam > ${sample}_head.sam
cat ${sample}_unique.sam >> ${sample}_head.sam
samtools sort -@ 20 -o ${sample}.bam ${sample}_head.sam
samtools index ${sample}.bam
rm ${sample}_unique.sam
rm ${sample}_head.sam
stringtie -p 20 -G path/to/Wm82.gtf -e -B -o ${sample}.gtf -A ${sample}.tsv ${sample}.bam
done

2.4 从tsv文件提取fpkm信息

ls *tsv|cut -d '.' -f 1|while read id;
do
nohup cat ${id}.tsv | awk -F "\t" '{print $1"\t"$8}' > fpkm/${id}.fpkm &
done

2.5 从gtf文件提取count信息

ls path/to/gtf/*gtf > gtf
ls path/to/fastq/*R2* | cut -d '/' -f 7 | cut -d '_' -f 1-3 > name
paste name gtf > prepDE.list
python path/to/prepDE.py -i path/to/prepDE.list -g prepDE_counts.csv -t prepDE_transcript.csv
sed 's/,/\t/' prepDE_counts.csv > all_count.txt

3、差异基因分析

DESeq2以基因count为输入文件

以15天的大豆根和根瘤为例：

3.1 DESeq2差异分析

将基因count文件下载到本地，使用R包DESeq2进行差异分析

rm(list=ls())
options(stringsAsFactors = F)
library(DESeq2)
library(tidyverse)

counts <- read_csv('path/to/D15_counts_prep.csv')
exprSet <- counts[,2:7]
rownames(exprSet) <- counts$gene

group_list <- factor(c('R15','R15','R15','N15','N15','N15'),
                     levels=c('R15','N15'))
colData <- data.frame(row.names=colnames(exprSet),
                      group_list=group_list)

dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds)
sizeFactors(dds)
res <- results(dds)
res_d <- subset(res,abs(log2FoldChange) > 1 & padj < 0.05)

DEG <- as.data.frame(res_d)
DEG <- na.omit(DEG)

D15_up <- DEG %>%
  filter(log2FoldChange>0)
D15_up <- tibble(rownames(D15_up))
names(D15_up) <- c("ID")
write_tsv(D15up, file='D15_up.txt')

D15_down <- DEG %>%
  filter(log2FoldChange<0)
D15_down <- tibble(rownames(D15_down))
names(D15_down) <- c("ID")
write_tsv(D15down, file='D15_down.txt')

#火山图
threshold <- as.factor(ifelse(DEG$padj < 0.05 & abs(DEG$log2FoldChange) >= 2,
                      ifelse(DEG$log2FoldChange >= 2 ,'Up','Down'),'Not'))
ggplot(DEG,aes(x=log2FoldChange,y=-log10(padj),colour=threshold))+
  geom_point(size=0.5,alpha=1)+
  labs(x='log2FC',y='-log10(qvalue)')+
  ylim(0,300) + xlim(-15,25)+
  scale_color_manual(values=c("#30A9DE",'grey',"#ff4e50"))+
  geom_vline(xintercept = c(-2, 2), lty = 2,colour="#000000")+
  geom_hline(yintercept = c(1), lty = 2,colour="#000000")+
  theme(panel.grid = element_blank(),
        axis.text=element_text(size=14),
        axis.title=element_text(size=14),
        legend.position = "none"
        )