数据的统计分布

通过定量PCR测量的不同细胞之间同一基因的表单水平服从 对数正太分布

使用负二项分布作为它们对RNA-seq计数的建模基础 然而 很大部分值为零，导致很难用该分布拟合（白说了）

最新的论文认为可以用 Poisson-Twdddie family 建模  对应的R包 tweeDESeq

重点工具：

推荐：DESeq 2

非参数方法 SAMSeq ，NOISeq

limma

归一化

大多软件在内部执行归一化

关于选择差异表达分析的软件包

决策树如下

选择要检验差异表达的特征的类型，用于：

差异表达的外显子：DEXSeq

差异表达的异构体：BitSeq,Cuffdiff 或 ebSeq

差异表达的基因：选择实验设计的类型

复杂的设计（多因素）：DESeq  edgeR  limma

组间的简单比较：多少个生物学重复？

每个组超过5个生物学重复：SAMSeq

少于5个：DESeq   edgeR   limma

 

典型工作流程

1 fastq -> 质量控制和预处理  ->

2 比对到基因组 或 转录组

2.1  比对到基因组 -> BAM  -> 读段计数 HTseq 、featureCounts 、 BEDTools等  -> 差异表达分析 DESeq edgeR  limma SAMSeq  DEXSeq 等

2.2  比对到 转录组 -> BAM -> 基于异构体反褶积的差异表达分析    Cuffdiff 、BitSeq 、ebSeq等

 

 

火山图

做火山图使用差异倍数 ford change 及 p值   代码如下

def test():
    file_url = "e://生物信息//测试组学111-ABC transportor.xlsx"
    df = DataFrame(pd.read_excel(file_url , sheet_name="2",index_col="Gene_ID"))

    #预处理
    #data.mean(axis = 1)
    df.eval("""
        AD1_avg = (AD1_1 + AD1_2 + AD1_3 )/3
        AD3_avg = (AD3_1 + AD3_2 + AD3_3 )/3
        AD5_avg = (AD5_1 + AD5_2 + AD5_3 )/3
        CK1_avg = (CK1_1 + CK1_2 + CK1_3 )/3
        CK3_avg = (CK3_1 + CK3_2 + CK3_3 )/3
        CK5_avg = (CK5_1 + CK5_2 + CK5_3 )/3
        AD_avg = (AD1_avg + AD3_avg + AD5_avg) / 3
        CK_avg = (CK1_avg + CK3_avg + CK5_avg) / 3
        FoldChange = AD_avg / CK_avg
        """, inplace=True)

    df.to_excel("e://生物信息//Volcano.xlsx" , sheet_name="1")

    df2 = DataFrame(df,columns=["FoldChange"])

    # fold change 通常为了防止取log2时产生NA，给表达值加1 图像中心向右移动到1
    df2 = np.log2(df2+1)

    # 和对照组 做t检验获得p值
    a = df[['AD1_avg','AD3_avg','AD5_avg']]
    b = df[['CK1_avg','CK3_avg','CK5_avg']]
    t,p = stats.ttest_rel(a, b, axis=1) #按行配对

    df2.insert(loc=1,column='pvalue',value=p)
    df2['log(pvalue)'] = -np.log10(df2['pvalue'])

    #选定的差异基因标准是 I.差异倍数的绝对值大于1，II. 差异分析的P值小于0.05
    df2['sig'] = 'normal'
    df2['size'] = np.abs(df2['FoldChange']) / 10
    df2.loc[(df2.FoldChange > 1) & (df2.pvalue < 0.05), 'sig'] = 'up'
    df2.loc[(df2.FoldChange < -1) & (df2.pvalue < 0.05), 'sig'] = 'down'
    print(df2)
    df2.to_excel("e://生物信息//Volcano.xlsx", sheet_name="火山图")
    ax = sns.scatterplot(x="FoldChange", y="log(pvalue)",
                         hue='sig',
                         hue_order=('down', 'normal', 'up'),
                         palette=("#377EB8", "grey", "#E41A1C"),
                         data=df2)
    ax.set_ylabel('-log(pvalue)', fontweight='bold')
    ax.set_xlabel('FoldChange', fontweight='bold')
    plt.show()

    #筛选差异基因
    fold = df2['FoldChange'].tolist()
    pvalue = df2['pvalue'].tolist()
    fold_cutoff = 1 #阈值
    pvalue_cutoff = 0.05 #显著性

    filtered_ids = []
    for i in range(0, len(fold)):
        if (abs(fold[i]) >= fold_cutoff) and (pvalue[i] <= pvalue_cutoff):
            filtered_ids.append(i)

    filtered = df2.iloc[filtered_ids, :]
    print(filtered)
    print("Number of DE genes: ")
    print(len(filtered.index))