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

image.png

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,743评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,296评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,285评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,485评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,581评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,821评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,960评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,719评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,186评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,516评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,650评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,329评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,936评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,757评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,991评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,370评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,527评论 2 349