搜索
查看: 2624|回复: 0

[mRNA-seq] 转录组作业(七)差异表达分析

[复制链接]

11

主题

22

帖子

330

积分

中级会员

Rank: 3Rank: 3

积分
330
发表于 2017-11-6 21:16:40 | 显示全部楼层 |阅读模式
导入数据







> setwd("D:/RNA_seq/aligned/matrix")
> options(stringsAsFactors = FALSE) #这一步很重要,一定要写FALSE
> control<- read.table("SRR3589956.count",sep="\t",col.names = c("gene_id","control"))
> View(control)
> rep1<-read.table("SRR3589957.count",sep="\t",col.names = c("gene_id","rep1"))
> rep2<-read.table("SRR3589958.count",sep="\t",col.names = c("gene_id","rep1"))
> raw_count <- merge(merge(control,rep1,by="gene_id"),rep2,by="gene_id")
> head(raw_count)
                 gene_id  control   rep1.x   rep1.y
1 __alignment_not_unique 11246816 11094117 11713427
2            __ambiguous  3486284  3824345  3296626
3           __no_feature  2233241  1740948  1418675
4          __not_aligned  1143671  1307561  1242118
5        __too_low_aQual   754973   738068   656393
6   ENSG00000000003.14_2     1650     1971     1686
> raw_count_filt <- raw_count[-1:-5,] #去除没有匹配的reads
> View(raw_count_filt)
> ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
这一步可以看一下得到的ENSENBL长什么样子
> row.names(raw_count_filt) <- ENSEMBL #把矩阵的name变为可查询的id
> View(raw_count_filt)
> summary(raw_count_filt)
   gene_id             control             rep1.x        
Length:60492       Min.   :     0.0   Min.   :     0.0  
Class :character   1st Qu.:     0.0   1st Qu.:     0.0  
Mode  :character   Median :     0.0   Median :     0.0  
                    Mean   :   728.6   Mean   :   757.2  
                    3rd Qu.:    29.0   3rd Qu.:    30.0  
                    Max.   :328729.0   Max.   :244930.0  
     rep1.y        
Min.   :     0.0  
1st Qu.:     0.0  
Median :     0.0  
Mean   :   648.1  
3rd Qu.:    20.0  
Max.   :205018.0  
> GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",] #提取数据
> GAPDH
                             gene_id control rep1.x rep1.y
ENSG00000111640 ENSG00000111640.14_2   99138 125066 128452
> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
                            gene_id control rep1.x rep1.y
ENSG00000105127 ENSG00000105127.8_2    2364   1074   1015

#数据填坑
>colnames(raw_count_filt)
[1] "gene_id" "control" "rep1.x"  "rep1.y"
> delta_mean <- abs(mean(raw_count_filt$rep1.x) - mean(raw_count_filt$rep1.y))
> sampleNum <- length(raw_count_filt$control)
> sampleMean <- mean(raw_count_filt$control)
> control2 <- integer(sampleNum)
> for (i in 1:sampleNum){
+ if(raw_count_filt$control < sampleMean){
+ control2 <- raw_count_filt$control + abs(raw_count_filt$rep1.x - raw_count_filt$rep1.y)
+ } else{
+ control2 <- raw_count_filt$control + rpois(1,delta_mean)
+ }
+ }





统计学知识
样品服从正态分布,两个样本可以用t检验,多个样本的比较可以用方差分析。芯片数据的资料往往认为是服从正态分布。差异分析就是把t检验或是方差分析应用到样品的每一个基因上面去。线性回归模型指的是自变量X和因变量Y为线性关系,广义线性模型 包括一般线性回归,贝叶斯回归,方差分析等。是一般线性模型的推广,首先自变量可以是离散的,也可以是连续的。离散的可以是0-1变量,也可以是多种取值的变量。与线性回归模型相比较,随机误差项不一定服从正态分布,可以服从二项,泊松,负二项,高斯分布等等
高通量测序得到的read count之间的差异服从负二项分布,负二向分布有两个参数,均值(mean)和离散值(dispersion). 离散值描述方差偏离均值的程度。泊松分布可以认为是负二向分布的离散值为1,也就是均值等于方差(mean=variance)的情况。



library(DESeq2)
countData<-raw_count_filt[,2:5]
condition<-factor(c("control","KD","KD","control"))
dds<-DESeqDataSetFromMatrix(countData,DataFrame(condition),design=~condition)
head(dds)
str(dds)
?DESeqDataSetFromMatrix
nrow(dds)
ncol(dds)
dds<-dds[rowSums(counts(dds))>1,]
nrow(dds)
dds<-DESeq(dds)
res<-results()
res<-results(dds)
mcols(res,use.names=TRUE)
summary(res)
# out of 31282 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 3294, 11%
LFC < 0 (down)   : 2453, 7.8%
outliers [1]     : 0, 0%
low counts [2]   : 14556, 47%
(mean count < 15)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?resultts

#MA图
plotMA(res,ylim=c(-5,5))
topGene<- rownames(res)[which.min(res$padj)]
with(res[topGene,],{
points(baseMean,log2FoldChange,col="dodgerblue",cex=2,lwd=2)
  text(baseMean,log2FoldChange,topGene,pos=2,col="dodgerblue")
})
#lfcShrink收缩
res.shrink<-lfcShrink(dds,contrast = c("condition","KD","control"),res=res)
plotMA(res.shrink,ylim=c(-5,5))
topGene<- rownames(res)[which.min(res$padj)]
with(res[topGene,],{
  points(baseMean,log2FoldChange,col="dodgerblue",cex=2,lwd=2)
  text(baseMean,log2FoldChange,topGene,pos=2,col="dodgerblue")
})

#找出差异基因
res.deseq2<-subset(res,padj<0.05)
source("http://Bioconductor.org/biocLite.R")
library(xlsx)
write.xlsx2(res.deseq2,"DEG.xlsx",colnames=T,row.names=T,append=T)












走一遍流程,代码以及全部知识都是来自徐洲更大神在生信媛的发文。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:Python小测003
下一篇:Python小测004
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-22 03:10 , Processed in 0.039813 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.