hisat2+stringtie+deseq2 RNAseq（二）

测序得到的是几百bp的短read，相当于把拼图打散了给你。如果没有参考基因组，从头组装会有很大的工作量，而且对于有参考基因组的动植物来说，不依靠参考基因组进行RNAseq 所得到的准确性较低，因此我们需要下载参考基因组，将测序得到的短序列与参考基因组进行比对。

因为我们是做的猪细胞测序，因此我们需要下载猪的参考基因组
首先是进入NCBI
下载猪的参考基因组

image.png

#下载参考基因组
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.fna.gz

#然而参考基因组是一部无字天书，要想解读书中的内容，需要额外的注释信息协助。

#下载注释文件
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.gff.gz

#注释文件一般以gff 或者gtf结尾

那么接下来就是将测序数据与参考基因组基因比对

cd /8Tdata/liuhailong/data/
for i in {1..18};
do
    hisat2 -p 8 --dta -x /8Tdata/liuhailong/genome/sus -1 R${i}_1.out.fq -2  R${i}_2.out.fq --rna-strandness RF -S R${i}.sam --known-splicesite-infile /8Tdata/liuhailong/splicesites.txt --novel-splicesite-outfile /8Tdata/liuhailong/splicesite.bed
done
#这里用的是循环比对，
-p 进程数
--dta 对齐比对
-x 参考基因组
-1 第一个测序数据
-2 第二个测序数据
-S 输出结果
--rna-standness 表示链特异性，如果没有就不管

将生成的sam 转化为bam
samtools sort -@ 2 -o SRR${i}.bam SRR${i}.sam

使用tringtie拼接序列
for i in {1..18};
do
    #samtools sort -@ 8 -o R${i}.bam R${i}.sam
    stringtie  -p 8  --rf -G /8Tdata/liuhailong/genome/Sus_scrofa.Sscrofa11.1.97.gtf  -o /8Tdata/liuhailong/final_data/ballgown/R${i}/R${i}.gtf -B -e  R${i}.bam

done

这个时候就要利用一个python文件
stringtie官方提供了一个python文件

image.png

利用prepDE.py就可以将生成的gtf文件转化为计数文件
生成DEseq2能够读取的read count 矩阵

python ~/Software/prepDE.py -i gtflist.txt -g countRes/gene_count.csv -t countRes/transcript.csv

附：gtflist.txt格式：

SRR3469478    ./SRR3469478.gtf

SRR3469479      ./SRR3469479.gtf

SRR4421540      ./SRR4421540.gtf

SRR4421541      ./SRR4421541.gtf
或者直接硬转化
python prepDE.py  -g /8Tdata/liuhailong/final_data/count/R${i}.gene_count_matrix.csv -t /8Tdata/liuhailong/final_data/ballgown/R${i}.transcript_count_matrix.csv
也是可以的

转化的结果如图所示，这样就可以直接用deseq2进行比对

image.png

library(DESeq2)

c1<-read.table("PK-1_L02_503.count",row.names=1,header=T)
c2<-read.table("PK-2_L02_504.count",row.names=1,header=T)
c3<-read.table("PK-3_L02_505.count",row.names=1,header=T)

c4<-read.table("PK+9_L04_509.count",row.names=1,header=T)
c5<-read.table("PK+10_L04_510.count",row.names=1,header=T)
c6<-read.table("PK+11_L03_511.count",row.names=1,header=T)
#将文件读入进去
counts<-cbind(c1,c2,c3,c4,c5,c6)
#合并文件
colnames(counts)<-c("503","504","505","509","510","511")
condition<-factor(c('P150','P150','P150','NC','NC','NC'))
#设置变量因子
coldata<-data.frame(row.names=colnames(counts),condition)

## 构建dds矩阵
dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design= ~condition)

#对原始dds进行normalize
dds2 <- DESeq(dds)
res <- results(dds2)

diff_gene_deseq2 <-subset(res,padj < 0.001 & (log2FoldChange > 2 | log2FoldChange < -2))
diff_gene_deseq2_out <- cbind(row.names(diff_gene_deseq2),diff_gene_deseq2$baseMean,diff_gene_deseq2$log2FoldChange,
diff_gene_deseq2$lfcSE,diff_gene_deseq2$stat,diff_gene_deseq2$pvalue,diff_gene_deseq2$padj)

colnames(diff_gene_deseq2_out)<-c('gene',colnames(diff_gene_deseq2))
resdata <-  merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
# 得到csv格式的差异表达分析结果
List<-NULL;
for(i in 1:dim(diff_gene_deseq2_out)[1]){
  gene<-diff_gene_deseq2_out[i,1]
  p<-resdata[which(resdata$Row.names==gene),]
  List<-rbind(List,p)
}
write.csv(List,file= "R123_vs_JEV.difference-gene.p0.001-f2.result.csv",row.names = F)

write.csv(resdata,file= "PK-MOCK_vs_PK+TGEV.all.result.csv",row.names = F)
这样即可得到需要的结果

image.png

接下来就是需要将富集的基因进行功能注释
在这里我们就推荐 clusterProfiler 对基因ID的转化

source("http://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
library(clusterProfiler)
library(org.Ss.eg.db)
keytypes(org.Mm.eg.db)   #用来看所能选择ID的键值
data<-read.table("GX_VS_HB.all.result.csv",header = T,sep = ",")
colnames(data)[1]<-c("ensembl_gene_id")   #修改Ensemble基因所在的位置
gene<-data$ensembl_gene_id  #提取基因
eg <- bitr(gene, fromType="ENSEMBL", toType=c("SYMBOL"), OrgDb="org.Mm.eg.db")
colnames(eg)[1]<-c("ensembl_gene_id")
diff_name<-merge(x=data,y=eg,by="ensembl_gene_id",all.x=TRUE) #合并
write.csv(diff_name,file= "GX_VS_HB.all.csv",row.names = F) #输出

go <- enrichGO(symbol, OrgDb = org.Bt.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'SYMBOL')
write.csv(go,file= "333.csv",row.names = F) #输出

kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'bta',
               pvalueCutoff = 0.05)

kk <- enrichKEGG(gene = gene,
                 organism ='hsa',
                 pvalueCutoff = 0.05,
                 qvalueCutoff = 0.1,
                 minGSSize = 1,
                 #readable = TRUE ,
                 use_internal_data =FALSE)

最后将富集的结果可视化即可

image.png

如果有什么问题可以联系作者