美文网首页
Deseq2 分析RNA-seq的reads counts矩阵

Deseq2 分析RNA-seq的reads counts矩阵

作者: 余绕 | 来源:发表于2020-06-18 20:13 被阅读0次

本数据为自己随意产生的,故分析时有提示!

library(DESeq2)
> rm(list=ls())
> count_data=read.table('DEseq.txt',header = TRUE,sep="\t")
> head(count_data)
  Gene_id ev_1 ev_2 sample_1 sample_2 sample_3
1   gene1   53   32       97       59       27
2   gene2   24    3       74        2       73
3   gene3   46    5       88       34       42
4   gene4   75   58       34       21        1
5   gene5   90   65       25        1       53
6   gene6   89   99       96       12       86
> row.names(count_data)=count_data$Gene_id#把第一列变成行名
> head(count_data)
      Gene_id ev_1 ev_2 sample_1 sample_2 sample_3
gene1   gene1   53   32       97       59       27
gene2   gene2   24    3       74        2       73
gene3   gene3   46    5       88       34       42
gene4   gene4   75   58       34       21        1
gene5   gene5   90   65       25        1       53
gene6   gene6   89   99       96       12       86
> count_tab=count_data[,-1] #去除第一列
> head(count_tab)
      ev_1 ev_2 sample_1 sample_2 sample_3
gene1   53   32       97       59       27
gene2   24    3       74        2       73
gene3   46    5       88       34       42
gene4   75   58       34       21        1
gene5   90   65       25        1       53
gene6   89   99       96       12       86
> sample_Data=read.table('Phynotype.txt',header = TRUE,sep="\t")
> sample_Data
  SampleID condition
1     ev_1        EV
2     ev_2        EV
3 sample_1    Sample
4 sample_2    Sample
5 sample_3    Sample
> sample_Data$condition=factor(sample_Data$condition,c("EV","Sample")) #设置分组顺序,这里因为E在S前,不设置也可以
> sample_Data
  SampleID condition
1     ev_1        EV
2     ev_2        EV
3 sample_1    Sample
4 sample_2    Sample
5 sample_3    Sample
> dds<- DESeqDataSetFromMatrix(countData = count_tab,
+                              colData =sample_Data,
+                              design = ~condition)
converting counts to integer mode
> dds<-DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
> sizefactor<-sizeFactors(dds) #获得各组的size factors
> sizefactor
    ev_1     ev_2 sample_1 sample_2 sample_3 
1.151989 1.185014 1.117790 1.141435 1.141129 
> resultsNames(dds)#list the coefficients,获得比较组
[1] "Intercept"              "condition_Sample_vs_EV"
> res<-results(dds,name="condition_Sample_vs_EV")
> head(res)
log2 fold change (MLE): condition Sample vs EV 
Wald test p-value: condition Sample vs EV 
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat            pvalue              padj
             <numeric>          <numeric>         <numeric>          <numeric>         <numeric>         <numeric>
gene1 47.0279585924932  0.566053769686641  0.76867337531233  0.736403507480197 0.461485170730465 0.993661998703688
gene2  31.058228262017   1.91337507463314  1.31773242322843   1.45202094211615 0.146495782309875 0.993661998703688
gene3 37.8939636752357   1.13461160222758  1.03127890545406   1.10019859441227 0.271245602893762 0.993661998703688
gene4 32.7481554501855  -1.78404566313186   1.1589998211077  -1.53929761734283  0.12373165537549 0.993661998703688
gene5 40.5328677696536  -1.51709659185057   1.0826512214148  -1.40127915790649 0.161130611211371 0.993661998703688
gene6 66.5123710070533 -0.489985663564088 0.995504546061613 -0.492198318433156 0.622579153093233 0.993661998703688
> res<-res[order(res$padj),] #按照padjust value排序
> head(res)
log2 fold change (MLE): condition Sample vs EV 
Wald test p-value: condition Sample vs EV 
DataFrame with 6 rows and 6 columns
                baseMean   log2FoldChange             lfcSE             stat               pvalue
               <numeric>        <numeric>         <numeric>        <numeric>            <numeric>
gene372 44.0121189310388 4.18981724397574 0.914268181788086 4.58270048923881 4.59009265192878e-06
gene206 28.0027885585295 5.75134738523254  1.39157480607446 4.13297751592436 3.58093768694533e-05
gene281 39.1423657824902 3.72547233692155 0.913615632247554 4.07772394147487 4.54787187421144e-05
gene490 46.5639545775845 3.83712142893487 0.940257744471528 4.08092509899137 4.48568095546619e-05
gene600 32.4274383805788 4.94359656208485  1.29720518270051 3.81095961380088 0.000138428383468241
gene504 29.2387021124829 4.78408878059009  1.36221357333734 3.51199611737048 0.000444754496479191
                       padj
                  <numeric>
gene372  0.0032543756902175
gene206 0.00806110289703977
gene281 0.00806110289703977
gene490 0.00806110289703977
gene600  0.0196291447757966
gene504  0.0525551563339577
> resDF<-as.data.frame(res)
> resDF$Gene_id=row.names(resDF)#创建geneID一列,产生的数据在最后一列
> head(resDF)
        baseMean log2FoldChange     lfcSE     stat       pvalue        padj Gene_id
gene372 44.01212       4.189817 0.9142682 4.582700 4.590093e-06 0.003254376 gene372
gene206 28.00279       5.751347 1.3915748 4.132978 3.580938e-05 0.008061103 gene206
gene281 39.14237       3.725472 0.9136156 4.077724 4.547872e-05 0.008061103 gene281
gene490 46.56395       3.837121 0.9402577 4.080925 4.485681e-05 0.008061103 gene490
gene600 32.42744       4.943597 1.2972052 3.810960 1.384284e-04 0.019629145 gene600
gene504 29.23870       4.784089 1.3622136 3.511996 4.447545e-04 0.052555156 gene504
> resDF<-resDF[,c(7,1:6)]#把最后一列放在最前面
> head(resDF)
        Gene_id baseMean log2FoldChange     lfcSE     stat       pvalue        padj
gene372 gene372 44.01212       4.189817 0.9142682 4.582700 4.590093e-06 0.003254376
gene206 gene206 28.00279       5.751347 1.3915748 4.132978 3.580938e-05 0.008061103
gene281 gene281 39.14237       3.725472 0.9136156 4.077724 4.547872e-05 0.008061103
gene490 gene490 46.56395       3.837121 0.9402577 4.080925 4.485681e-05 0.008061103
gene600 gene600 32.42744       4.943597 1.2972052 3.810960 1.384284e-04 0.019629145
gene504 gene504 29.23870       4.784089 1.3622136 3.511996 4.447545e-04 0.052555156
> final_data=merge(resDF,count_data,by="Gene_id") #合并数据,把处理前的readscounts合并进去
> head(final_data)
  Gene_id baseMean log2FoldChange     lfcSE       stat    pvalue     padj ev_1 ev_2 sample_1 sample_2 sample_3
1   gene1 47.02796      0.5660538 0.7686734  0.7364035 0.4614852 0.993662   53   32       97       59       27
2  gene10 29.79123      0.6988787 1.2667329  0.5517175 0.5811419 0.993662   49    1       71       20       28
3 gene100 42.28342     -0.1866548 0.8440605 -0.2211391 0.8249841 0.993662   73   33       56       18       62
4 gene101 48.47256     -0.1330532 0.7450646 -0.1785794 0.8582679 0.993662   80   39       85       34       39
5 gene102 35.12121      0.7993091 0.9708976  0.8232682 0.4103555 0.993662   22   35       43       19       82
6 gene103 48.32332      0.7512332 0.6814179  1.1024560 0.2702635 0.993662   43   37       72       45       79
>

相关文章

网友评论

      本文标题:Deseq2 分析RNA-seq的reads counts矩阵

      本文链接:https://www.haomeiwen.com/subject/qfrpxktx.html