RNA-seq数据分析
5.43K
0
数据的统计分布
通过定量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))