02基因差异表达分析

R包安装

使用R3.6.0，进行

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ballgown")
##也可以直接使用library加载
library(ballgown)

报错

BiocManager在安装过程中没有权限

https://stackoverflow.com/questions/41839214/installation-path-not-writable-r-unable-to-update-packages 还是没有解决

# Installation path not writeable
installed.packages()[, c("Package", "LibPath")]
#查看这些包的安装路径，直接改变权限

1.加载表达文件

当stringtie生成的*.ctab文件都在同一个目录下的时候，直接使用样本所在的父目录作为参数。samplePattern样本的前缀，这个例子中包含有20个样本的数据

# make the ballgown object:
#根据目录来提取表达
bg = ballgown(dataDir=data_directory, samplePattern='sample', meas='all')
bg
#直接指定样本文件夹所在路径
ballgown(samples=c(文件目录),meas='all')
## ballgown instance with 100 transcripts and 20 samples

2.使用插槽访问里面的各种数据

*expr(ballgown_object_name, <EXPRESSION_MEASUREMENT>)
#例如
transcript_fpkm = texpr(bg, 'FPKM')
transcript_cov = texpr(bg, 'cov')
whole_tx_table = texpr(bg, 'all')
gene_expression=gexpr(bg) ##制作多个样本的表达谱

差异表达分析

Ballgown provides a wide selection of simple, fast statistical methods for testing whether transcripts are differentially expressed between experimental conditions or across a continuous covariate (such as time).

总共15个分类，每类里有两个重复

##指定每个样本的分组信息，总共15个样，每个样有两个重复
pData(bg) = data.frame(id=sampleNames(bg), group=rep(seq(0,14), each=2))
##进行两两间的差异表达分析,对第一组和第二组进行差异表达分析
##这里将所有的样品随机分成两类提，进行差异表达分析
sex = sample(c('M','F'), size=nrow(pData(bg)), replace=TRUE)
age = sample(21:52, size=nrow(pData(bg)), replace=TRUE)

mod = model.matrix(~ sex + age + pData(bg)$group + pData(bg)$time)
mod0 = model.matrix(~ pData(bg)$group + pData(bg)$time)
adjusted_results = stattest(bg, feature='transcript', meas='FPKM', mod0=mod0, mod=mod)

也可以得到matrix用于到其他模型中,例如DEseq2等。

分析两个重复间的重复性和差异表达基因

##使用ballgown读取表达量
bg=ballgown(samples=c("1-1-1","1-1-2"),,meas='all')
##合并成一个表达谱
gene_expression=gexpr(bg)
##计算相关性
cor.test(gene_expression[,1],gene_expression[,2])
##计算重复间差异表达的基因
pData(bg)=data.frame(id=sampleNames(bg),group=rep(seq(0,1),each=2))
stat_results=stattest(bg, feature='transcript', meas='FPKM', covariate='group')

使用DESeq2进行差异表达分析

htseq统计比对到gene的read数目

for k in ${all_sample[@]};
do
bsub -q smp -n 10 -R span[hosts=1] -J ${k} -e ${k}.err -o ${k}.out "
 htseq-count  -f bam -n 10    -r pos -i gene_id  -m union -q  ${inputDir}/${k}/${k}_sort.bam   /data/cotton/zhenpingliu/genome/genome_data/Ghirsutum_genome_HAU_v1.1/Ghirsutum_gene_model.gtf   >${k}_count.txt"
sleep 1
done

进行差异表达分析

gene_expression=as.matrix(filter_gene[,-1])
row.names(gene_expression)=filter_gene[,1]
cond <- factor(rep(1:2, each=2))
dds=DESeqDataSetFromMatrix(gene_expression,DataFrame(cond),~cond)
# 指定哪一组作为control
dds$cond <- relevel(dds$cond, ref = "2")
##进行差异表达分析
dds <- DESeq(dds)
res <- results(dds)
##输出到文件
res$geneId=rownames(res)
res=res[c('geneId','baseMean','log2FoldChange','lfcSE','stat','pvalue','padj')]
write.table(res,file = "test",quote = FALSE,sep="\t",row.names = F)

##统计差异表达基因的数目
awk '$3>=1||$3<=-1{print $0}' difference_expression.txt |awk '$6<=0.05{print $0}'|awk '{print $1"\t"$3}'|wc -l

绘制火山图

##将DESeq2数据转为data.frame
res=as.data.frame(res)
##对上调、下调和没有变化的基因进行分类
res=  mutate(res,group=pvalue)
resDown=subset(res,log2FoldChange<=-1&pvalue<=0.05)
resDown$group='down'
resUp=subset(res,log2FoldChange>=1&pvalue<=0.05)
resUp$group='up'
resnoDiffer=subset(res,pvalue>0.05)
resnoDiffer$group='noDiff'
res2=rbind(resDown,resUp,resnoDiffer)
##绘制火山图和指明突变基因的位置

火山图

比较分类

• 1、3 和J在不同发育阶段比较

• 6、7和J在不同发育阶段比较

• 1、2和6、7在不同发育阶段比较

• 1、3在三个时期进行轨迹的比较

• 6、7在三个时期进行轨迹的比较

• J在三个时间进行轨迹的比较

参考

1. 转录组差异表达分析 https://github.com/alyssafrazee/ballgown

2. 官方pdf https://bioconductor.org/packages/release/bioc/manuals/ballgown/man/ballgown.pdf