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

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

使用負二項分布作為它們對RNA-seq計數(shù)的建?；A 然而 很大部分值為零，導致很難用該分布擬合（白說了）

最新的論文認為可以用 Poisson-Twdddie family 建模? 對應的R包 tweeDESeq

重點工具：

推薦：DESeq 2

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

limma

歸一化

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

關于選擇差異表達分析的軟件包

決策樹如下

選擇要檢驗差異表達的特征的類型，用于：

差異表達的外顯子：DEXSeq

差異表達的異構體：BitSeq,Cuffdiff 或 ebSeq

差異表達的基因：選擇實驗設計的類型

復雜的設計（多因素）：DESeq? edgeR? limma

組間的簡單比較：多少個生物學重復？

每個組超過5個生物學重復：SAMSeq

少于5個：DESeq? ?edgeR? ?limma

 

典型工作流程

1 fastq -> 質量控制和預處理? ->

2 比對到基因組 或 轉錄組

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

2.2? 比對到 轉錄組 -> BAM -> 基于異構體反褶積的差異表達分析? ? Cuffdiff 、BitSeq 、ebSeq等

 

 

火山圖

做火山圖使用差異倍數(shù) 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.差異倍數(shù)的絕對值大于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))