數(shù)據(jù)的統(tǒng)計(jì)分布

通過(guò)定量PCR測(cè)量的不同細(xì)胞之間同一基因的表單水平服從 對(duì)數(shù)正太分布

使用負(fù)二項(xiàng)分布作為它們對(duì)RNA-seq計(jì)數(shù)的建模基礎(chǔ) 然而 很大部分值為零，導(dǎo)致很難用該分布擬合（白說(shuō)了）

最新的論文認(rèn)為可以用 Poisson-Twdddie family 建模? 對(duì)應(yīng)的R包 tweeDESeq

重點(diǎn)工具：

推薦：DESeq 2

非參數(shù)方法 SAMSeq ，NOISeq

limma

歸一化

大多軟件在內(nèi)部執(zhí)行歸一化

關(guān)于選擇差異表達(dá)分析的軟件包

決策樹(shù)如下

選擇要檢驗(yàn)差異表達(dá)的特征的類(lèi)型，用于：

差異表達(dá)的外顯子：DEXSeq

差異表達(dá)的異構(gòu)體：BitSeq,Cuffdiff 或 ebSeq

差異表達(dá)的基因：選擇實(shí)驗(yàn)設(shè)計(jì)的類(lèi)型

復(fù)雜的設(shè)計(jì)（多因素）：DESeq? edgeR? limma

組間的簡(jiǎn)單比較：多少個(gè)生物學(xué)重復(fù)？

每個(gè)組超過(guò)5個(gè)生物學(xué)重復(fù)：SAMSeq

少于5個(gè)：DESeq? ?edgeR? ?limma

 

典型工作流程

1 fastq -> 質(zhì)量控制和預(yù)處理? ->

2 比對(duì)到基因組 或 轉(zhuǎn)錄組

2.1? 比對(duì)到基因組 -> BAM? -> 讀段計(jì)數(shù) HTseq 、featureCounts 、 BEDTools等? -> 差異表達(dá)分析 DESeq edgeR? limma SAMSeq? DEXSeq 等

2.2? 比對(duì)到 轉(zhuǎn)錄組 -> BAM -> 基于異構(gòu)體反褶積的差異表達(dá)分析? ? Cuffdiff 、BitSeq 、ebSeq等

 

 

火山圖

做火山圖使用差異倍數(shù) ford change 及 p值? ?代碼如下

def test():
    file_url = "e://生物信息//測(cè)試組學(xué)111-ABC transportor.xlsx"
    df = DataFrame(pd.read_excel(file_url , sheet_name="2",index_col="Gene_ID"))

    #預(yù)處理
    #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時(shí)產(chǎn)生NA，給表達(dá)值加1 圖像中心向右移動(dòng)到1
    df2 = np.log2(df2+1)

    # 和對(duì)照組 做t檢驗(yàn)獲得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) #按行配對(duì)

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

    #選定的差異基因標(biāo)準(zhǔn)是?I.差異倍數(shù)的絕對(duì)值大于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))