葡萄的參考基因組下載自NCBI，下載鏈接是
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/745/GCF_000003745.3_12X/

基本流程是

• Hiast2 比對(duì)
• samtools sam 轉(zhuǎn)bam
• stringtie組裝轉(zhuǎn)錄本
• gffcompare將stringtie輸出的gtf文件與參考基因組的注釋文件做比較得到一- 個(gè)merged.combine.gtf
• 使用merged.combine.gtf 這個(gè)文件對(duì)每個(gè)樣本計(jì)算表達(dá)量，輸出文件存儲(chǔ)到ballgown文件夾下，這一步用到的命令是 stringtie -e -B -p 8 -G merged.combined.gtf -o ballgown/L01/L01.gtf output_bam/L01.sorted.bam

image.png

• 接下來是R語言的ballgown包讀入數(shù)據(jù)獲取基因和轉(zhuǎn)錄本的表達(dá)量
  代碼是

library(ballgown)
library(genefilter)
library(dplyr)

pheno_data<-read.csv("pheno_data.txt",header=T)
grape_expr <- ballgown(dataDir = "ballgown",
                    samplePattern = "L0",
                    pData = pheno_data)

image.png

image.png

這一步可以拿到gene_id還有g(shù)ene_name ,FPKM的表達(dá)量，cov對(duì)用的應(yīng)該是reads count吧。

接下來是差異表達(dá)分析

代碼是

grape_expr_filter<-subset(grape_expr,
                          "rowVars(texpr(grape_expr))>1",
                          genomesubset=T)
results_genes <- stattest(grape_expr_filter,
                          feature = "gene",
                          covariate = "time_point",
                          getFC = TRUE,
                          meas = "FPKM")
#results_genes <- arrange(results_genes,pval)
results_genes%>%
  filter(abs(fc)>=2&pval<0.05) -> results_genes_diff

dim(results_genes_diff)
head(results_genes_diff)

現(xiàn)在有了基因id

image.png

接下來是使用clusterProfiler做go注釋

這里參考
https://guangchuangyu.github.io/cn/2017/07/clusterprofiler-maize/#disqus_thread

首先把這個(gè)基因id對(duì)應(yīng)著轉(zhuǎn)換成 ENTREZID ，這里需要借助對(duì)應(yīng)的gtf注釋文件

這里只關(guān)注蛋白編碼基因

grep 'gene_biotype "protein_coding"' 12X_genomic.gtf > 12X_protein_coding.gtf
#python
from gtfparse import read_gtf
known_proteincoding = read_gtf("12X_protein_coding.gtf")
known_proteincoding.to_csv("all_protein_coding.csv")

GO富集分析的R語言代碼

require(AnnotationHub)
hub<-AnnotationHub() #這一步對(duì)網(wǎng)路有要求
# aa<-query(hub,'zea')
# aa$title
# query(hub,'Malus domestica')
bb<-query(hub,"Vitis vinifera")
#bb$title
grape<-hub[['AH85815']]
# length(keys(grape))
# columns(grape)
protein_coding_all<-read.csv("all_protein_coding.csv",header=T)
df<-merge(results_genes_diff,grape_expr_filter@expr$trans,by.x="id",by.y="gene_id")
df1<-merge(df,protein_coding_all,by.x="gene_name",by.y="gene_id")
dim(df1)
gene_ids<- 
  df1$db_xref[!duplicated(df1$db_xref)]
gene_ids<-stringr::str_replace(gene_ids,"GeneID:","")


library(clusterProfiler)
bitr(keys(grape)[2],'ENTREZID',c("REFSEQ","GO",
                                 "ONTOLOGY","GENENAME",
                                 "SYMBOL"),grape)
res = enrichGO(gene_ids, 
               OrgDb=grape, pvalueCutoff=0.05, qvalueCutoff=0.05)
help(package="clusterProfiler")
dotplot(res)

最后的結(jié)果是

image.png

歡迎大家關(guān)注我的公眾號(hào)

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡(jiǎn)單小例子；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！