【生物信息学】基因差异分析Deg(数据读取、数据处理、差异分析、结果可视化)

news2024/11/27 5:37:55

目录

一、实验介绍

二、实验环境

1. 配置虚拟环境

2. 库版本介绍

3. IDE

三、实验内容

0. 导入必要的工具包

1. 定义一些阈值和参数

2. 读取数据

normal_data.csv部分展示

tumor_data.csv部分展示

3. 绘制箱型图

4. 删除表达量低于阈值的基因

5. 计算差异显著的基因

6. 初始化结果数据表格

7. 进行秩和检验和差异倍数计算

8. 可视化分析

9. 代码整合


一、实验介绍

        本实验完成了基因差异分析,包括数据读取、数据处理( 绘制箱型图、删除表达量低于阈值的基因、计算差异显著的基因)、差异分析(进行秩和检验和差异倍数计算)等,成功识别出在正常样本与肿瘤样本之间显著表达差异的基因,并对其进行了进一步的可视化分析(箱型图、差异倍数fold分布图、热力图和散点图)。

        基因差异分析是研究不同条件下基因表达差异的重要手段,能够帮助我们理解生物体内基因调控的变化及其与表型特征的关联。本实验旨在探索正常样本与肿瘤样本之间基因表达的差异,并识别差异显著的基因。

二、实验环境

    本系列实验使用了PyTorch深度学习框架,相关操作如下(基于深度学习系列文章的环境):

1. 配置虚拟环境

深度学习系列文章的环境

conda create -n DL python=3.7 
conda activate DL
pip install torch==1.8.1+cu102 torchvision==0.9.1+cu102 torchaudio==0.8.1 -f https://download.pytorch.org/whl/torch_stable.html
conda install matplotlib
conda install scikit-learn

新增加

conda install pandas
conda install seaborn
conda install networkx
conda install statsmodels
pip install pyHSICLasso

注:本人的实验环境按照上述顺序安装各种库,若想尝试一起安装(天知道会不会出问题)

2. 库版本介绍

软件包本实验版本目前最新版
matplotlib3.5.33.8.0
numpy1.21.61.26.0
python3.7.16
scikit-learn0.22.11.3.0
torch1.8.1+cu1022.0.1
torchaudio0.8.12.0.2
torchvision0.9.1+cu1020.15.2

新增

networkx2.6.33.1
pandas1.2.32.1.1
pyHSICLasso1.4.21.4.2
seaborn0.12.20.13.0
statsmodels0.13.50.14.0

3. IDE

        建议使用Pycharm(其中,pyHSICLasso库在VScode出错,尚未找到解决办法……)

win11 安装 Anaconda(2022.10)+pycharm(2022.3/2023.1.4)+配置虚拟环境_QomolangmaH的博客-CSDN博客icon-default.png?t=N7T8https://blog.csdn.net/m0_63834988/article/details/128693741

三、实验内容

0. 导入必要的工具包

import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
from scipy.stats import ranksums
from pyHSICLasso import HSICLasso
from sklearn.preprocessing import LabelEncoder

1. 定义一些阈值和参数

p_cutoff = 0.001
FC_cutoff = 3
num_cutoff = 10
  • p_cutoff 是判断差异显著性的 p 值阈值。
  • FC_cutoff 是判断差异的倍数阈值。
  • num_cutoff 是基因表达量筛选的阈值。

2. 读取数据

ndata = pd.read_csv("normal_data.csv", index_col=0, header=0)
tdata = pd.read_csv("tumor_data.csv", index_col=0, header=0)

        从文件 "normal_data.csv" 和 "tumor_data.csv" 中读取正常样本和肿瘤样本的基因表达数据。

normal_data.csv部分展示

TCGA-BL-A13J-11A-13R-A10U-07TCGA-BT-A20N-11A-11R-A14Y-07TCGA-BT-A20Q-11A-11R-A14Y-07TCGA-BT-A20R-11A-11R-A16R-07TCGA-BT-A20U-11A-11R-A14Y-07TCGA-BT-A20W-11A-11R-A14Y-07TCGA-BT-A2LA-11A-11R-A18C-07TCGA-BT-A2LB-11A-11R-A18C-07TCGA-CU-A0YN-11A-11R-A10U-07TCGA-CU-A0YR-11A-13R-A10U-07TCGA-GC-A3BM-11A-11R-A22U-07TCGA-GC-A3WC-11A-11R-A22U-07TCGA-GC-A6I3-11A-11R-A31N-07TCGA-GD-A2C5-11A-11R-A180-07TCGA-GD-A3OP-11A-11R-A220-07TCGA-GD-A3OQ-11A-21R-A220-07TCGA-K4-A3WV-11A-21R-A22U-07TCGA-K4-A54R-11A-11R-A26T-07TCGA-K4-A5RI-11A-11R-A28M-07
?10000000000000000000
?22.55213.89335.459416.458334.42645.221911.34638.2564.30548.747121.00843.183706.012817.762140.84528.308813.72161.5099
?38.51195.730510.44177.68119.0621.35343.421925.75141.48346.54719.65152.64725.312110.733719.36419.13819.527318.536414.4093
?4117.4995103.9108148.8251100.9537162.7054128.9808118.3924124.8423180.5137131.3118127.6985145.632795.2855119.1427149.7972128.4473142.8571128.644139.1199
?5619.5833738.4077679.3286786.70361132.5581032.8771158.651143.321687.4096690.5882697.5129697.37611551.793556.22011018.907646.4851895.4832903.8232836.9674
?60000000000000000000
?7128.3422115.4856119.258120.6965120.930254.794582.700491.372966.5702153.5294259.3317272.8863291.500794.4976250.6016253.2622313.5504270.6093187.172
?80.7376000.39570.77520.547900.463800000000.3332000
?900000000000000.398700000

tumor_data.csv部分展示

TCGA-BT-A42F-01A-11R-A23W-07TCGA-C4-A0EZ-01A-21R-A24X-07TCGA-C4-A0F0-01A-12R-A10U-07TCGA-C4-A0F1-01A-11R-A034-07TCGA-C4-A0F6-01A-11R-A10U-07TCGA-C4-A0F7-01A-11R-A084-07TCGA-CF-A1HR-01A-11R-A13Y-07TCGA-CF-A1HS-01A-11R-A13Y-07TCGA-CF-A27C-01A-11R-A16R-07TCGA-CF-A3MF-01A-12R-A21D-07TCGA-CF-A3MG-01A-11R-A20F-07TCGA-CF-A3MH-01A-11R-A20F-07TCGA-CF-A3MI-01A-11R-A20F-07TCGA-CF-A47S-01A-11R-A23W-07TCGA-CF-A47T-01A-11R-A23W-07TCGA-CF-A47V-01A-11R-A23W-07TCGA-CF-A47W-01A-11R-A23W-07TCGA-CF-A47X-01A-31R-A23W-07TCGA-CF-A47Y-01A-11R-A23W-07TCGA-CF-A5U8-01A-11R-A28M-07
?100.367100000000002.1101000000.57370
?23.76137.328105.35982.80936.477713.93723.568440.37629.09149.08534.74093.9742.39047.021111.2489019.57226.34546.4641
?32.89595.52180.663604.5452.89226.81445.397541.06653.62386.080612.13348.68657.80335.667717.39535.188717.411218.324713.092
?4117.7889443.1501144.144692.6276136.3486238.55791.604211.2791100.480927.5497101.9194105.91878.537120.6983212.084670.1299106.1415150.7324103.2415149.5243
?51754.6361546.3974391.827868.8195685.42011752.1671229.951456.665966.06551691.1261220.8531279.2292167.0481399.5921939.577978.75961754.7171110.2251262.7651477.273
?600000000000000000000
?7192.1065306.195515.926866.9972157.381932.169973.471749.6115426.9924745.960348.8152402.5713229.9982266.8196120.8459472.8729135.8491640.3191553.0694511.0994
?802.9371000.7354000.59771.16351.059600.803500.50970.60420.485501.45034.58980.5285
?900000000000000000000

3. 绘制箱型图

ndata.iloc[1000:1010, ].transpose().plot(kind='box', title='Normal Sample Gene Boxplot', rot=30)
plt.tick_params(labelsize=10)
plt.savefig('Normal_Sample_Gene_box_plot', bbox_inches='tight')
plt.show()

        绘制正常样本中部分基因的箱型图,并保存为图片文件。

4. 删除表达量低于阈值的基因

ndata = ndata.iloc[29:, :]
tdata = tdata.iloc[29:, :]

        删除正常样本和肿瘤样本中表达量低于阈值的基因。

5. 计算差异显著的基因

gene_sig = ndata[ndata.mean(axis=1) > num_cutoff].index.intersection(tdata[tdata.mean(axis=1) > num_cutoff].index)
Ndata = ndata.loc[gene_sig]
Tdata = tdata.loc[gene_sig]

        选择正常样本和肿瘤样本中表达量高于阈值的基因,并保存为新的数据。

6. 初始化结果数据表格

p_value = [1.] * len(Ndata)
log2FoldChange = []
label = ['0'] * len(Ndata)
result = pd.DataFrame([p_value, log2FoldChange, label])
result = result.transpose()
result.columns = ['p_value', 'log2FC', 'label']
result.index = Ndata.index.tolist()
p = []

        创建一个结果数据表格,包含 p 值、log2 fold change 和标签,并初始化为默认值。

7. 进行秩和检验和差异倍数计算

for i in Ndata.index:
    p1 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='greater')[1]
    p2 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='less')[1]
    p3 = ranksums(Tdata.loc[i], Ndata.loc[i])[1]
    result.loc[i, 'log2FC'] = np.log2(np.average(Tdata.loc[i]) / np.average(Ndata.loc[i]))
    p.append(p3)
    if (p1 <= p_cutoff):
       result.loc[i, 'p_value'] = p1
       if result.loc[i, 'log2FC'] > np.log2(FC_cutoff):
           result.loc[i, 'label'] = '1'

    if (p2 <= p_cutoff):
       result.loc[i, 'p_value'] = p2
       if result.loc[i, 'log2FC'] < -np.log2(FC_cutoff):
           result.loc[i, 'label'] = '-1'

        对每个基因进行秩和检验,计算 p 值和差异倍数,并根据设定的阈值确定差异显著的基因,并给予标签。

8. 可视化分析

print('finished')
plt.hist(result['log2FC'], bins=10, color='blue', alpha=0.6, edgecolor='black')
plt.title("fold change")
plt.show()

9. 代码整合

import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
from scipy.stats import ranksums
from pyHSICLasso import HSICLasso
from sklearn.preprocessing import LabelEncoder


p_cutoff = 0.001
FC_cutoff = 3
num_cutoff = 10

# 读取数据ndata normal data,tdata tumor data
ndata = pd.read_csv("normal_data.csv", index_col=0, header=0)
tdata = pd.read_csv("tumor_data.csv", index_col=0, header=0)

# 箱型图
# 中位数
# QL称为下四分位数,表示全部观察值中有四分之一的数据取值比它小;
# QU称为上四分位数,表示全部观察值中有四分之一的数据取值比它大;
# IQR称为四分位数间距,是上四分位数QU与下四分位数QL之差,其间包含了全部观察值的一半。
# 异常值通常被定义为小于QL-1.5IQR或大于QU+1.5IQR的值。
ndata.iloc[1000:1010, ].transpose().plot(kind='box', title='Normal Sample Gene Boxplot', rot=30)
plt.tick_params(labelsize=10)
plt.savefig('Normal_Sample_Gene_box_plot', bbox_inches='tight')
plt.show()

# delete gene expression < num
ndata = ndata.iloc[29:, :]
tdata = tdata.iloc[29:, :]
gene_sig = ndata[ndata.mean(axis=1) > num_cutoff].index.intersection(tdata[tdata.mean(axis=1) > num_cutoff].index)
Ndata = ndata.loc[gene_sig]
Tdata = tdata.loc[gene_sig]
p_value = [1.] * len(Ndata)
# log2 fold change的意思是取log2,这样可以可以让差异特别大的和差异比较小的数值缩小之间的差距
log2FoldChange = []
label = ['0'] * len(Ndata)

result = pd.DataFrame([p_value, log2FoldChange, label])
result = result.transpose()
result.columns = ['p_value', 'log2FC', 'label']
result.index = Ndata.index.tolist()
p = []

# 秩和检验(也叫Mann-Whitney U-test),检验两组数据是否来自具有相同中位数的连续分布,检验它们是否有显著差异。
# ’less‘ 基于 x 的分布小于基于 y 的分布
# ‘greater’ x 的分布大于 y 的分布。

for i in Ndata.index:
    p1 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='greater')[1]
    p2 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='less')[1]
    p3 = ranksums(Tdata.loc[i], Ndata.loc[i])[1]
    result.loc[i, 'log2FC'] = np.log2(np.average(Tdata.loc[i]) / np.average(Ndata.loc[i]))
    p.append(p3)
    if (p1 <= p_cutoff):
       result.loc[i, 'p_value'] = p1
       if result.loc[i, 'log2FC'] > np.log2(FC_cutoff):
           result.loc[i, 'label'] = '1'

    if (p2 <= p_cutoff):
       result.loc[i, 'p_value'] = p2
       if result.loc[i, 'log2FC'] < -np.log2(FC_cutoff):
           result.loc[i, 'label'] = '-1'

print('finished')
# 查看差异倍数fold分布
plt.hist(result['log2FC'], bins=10, color='blue', alpha=0.6, edgecolor='black')
plt.title("fold change")
plt.show()


# result.to_csv('result.csv')
up_data_n = Ndata.loc[result.query("(label == '1')").index]
up_data_t = Tdata.loc[result.query("(label == '1')").index]
down_data_n = Ndata.loc[result.query("(label == '-1')").index]
down_data_t = Tdata.loc[result.query("(label == '-1')").index]
data = pd.concat([pd.concat([up_data_n, up_data_t], axis=1), pd.concat([down_data_n, down_data_t], axis=1)], axis=0)
deg_gene = pd.DataFrame(data.index.tolist())
deg_gene.to_csv('deg_gene.csv')
z1 = zscore(np.log2(data+1), axis=0)

sns.heatmap(data=z1, cmap="coolwarm", xticklabels=False)
plt.savefig('heatmap_plot', bbox_inches='tight')
plt.show()


p = -np.log10(p)
ax = sns.scatterplot(x="log2FC", y=p,
                     hue='label',
                     hue_order=('-1', '0', '1'),
                     palette=("#377EB8","grey","#E41A1C"),
                     data=result)
ax.set_ylabel('-log(pvalue)', fontweight='bold')
ax.set_xlabel('FoldChange', fontweight='bold')

plt.savefig('deg_p_fc_plot', bbox_inches='tight')
plt.show()


 

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1052988.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

Git多账号管理通过ssh 公钥的方式,git,gitlab,gitee

按照目前国内访问git&#xff0c;如果不科学上网&#xff0c;我们很大可能访问会超时。基于这个&#xff0c;所以我现在的git 配置已经增加到了3个了 一个公司gitlab&#xff0c;一个git&#xff0c;一个gitee. 以下基于这个环境&#xff0c;我们来说明下如何创建配置ssh公钥。…

VBA技术资料MF63:遍历形状并改变颜色

【分享成果&#xff0c;随喜正能量】人生&#xff0c;一站有一站的风景&#xff0c;一岁有一岁的味道&#xff0c;你的年龄应该成为你生命的勋章而不是你伤感的理由。生活嘛&#xff0c;慢慢来&#xff0c;你又不差&#xff01;。 我给VBA的定义&#xff1a;VBA是个人小型自动…

【数据结构与算法】栈与队列相关算法的实现

目录 检查括号是否成对出现 反转字符串 循环队列的实现 使用队列实现栈 使用栈实现队列 检查括号是否成对出现 算法要求 给定一个只包括 (&#xff0c;)&#xff0c;{&#xff0c;}&#xff0c;[&#xff0c;] 的字符串&#xff0c;判断该字符串是否有效。 有效字符串需满…

树莓派基本配置(2)

安装motion $sudo apt-get update $sudo apt-get install motion配置motion sudo nano /etc/default/motionsudo nano /etc/motion/motion.conf主要改这些参数 //让Motion作为守护进程运行 daemon on ... //用这个端口号来读取数据 stream_port 8081 ... //网络上其它主机…

2023年上海市安全员B证证模拟考试题库及上海市安全员B证理论考试试题

题库来源&#xff1a;安全生产模拟考试一点通公众号小程序 2023年上海市安全员B证证模拟考试题库及上海市安全员B证理论考试试题是由安全生产模拟考试一点通提供&#xff0c;上海市安全员B证证模拟考试题库是根据上海市安全员B证最新版教材&#xff0c;上海市安全员B证大纲整理…

正点原子嵌入式linux驱动开发——TF-A使用

上一篇笔记STM32MP157芯片的开发环境&#xff0c;之后就直接简写为MP1。为了保证安全ARM推出了 Arm Trusted Firmware的可信固件&#xff0c;简称 TF-A。它是一个开源的软件&#xff0c;最早是用在Armv8-A&#xff0c;ST也在MP1里面使用到了TF-A。它的作用就是隔离硬件&#xf…

递推+模拟---想好如何存储?

递推模拟 输入输出问题CCF-CSP考试历年真题题型分类分组输入——可能有多组测试数据&#xff0c;对于每组数据 递推---从前面已知态--->后续未知态AcWing 3777. 砖块AcWing 1208. 翻硬币AcWing 1211. 蚂蚁感冒AcWing 3433. 吃糖果AcWing 821. 跳台阶 模拟202212-2-csp-训练计…

C语言--atoi函数详解及模拟实现

本篇概要 本篇博客主要讲述atoi函数的定义&#xff0c;用法及模拟实现。 文章目录 本篇概要1.atoi函数简介2.atoi函数的原型及参数3.atoi函数的使用举例4.atoi函数的模拟实现 1.atoi函数简介 atoi() 是 C语言的一个标准库函数&#xff0c;定义在<stdlib.h>头文件中。 at…

【RocketMQ】基本使用:Java操作RocketMQ(rocketmq-client)

【RocketMQ】基本使用&#xff1a;Java操作RocketMQ&#xff08;rocketmq-client&#xff09; 1.引入依赖 <dependency><groupId>org.apache.rocketmq</groupId><artifactId>rocketmq-client</artifactId><version>4.3.2</version>…

公司知识库搭建步骤,知识库建设与运营的四个步骤分享

在知识管理方面&#xff0c;团队中的每一员&#xff0c;都像是一名独行侠&#xff0c;自己的知识&#xff0c;满足自己的需要&#xff0c;这其中&#xff0c;就造成了很多无意义的精力消耗。 公司知识库搭建必要性 比如&#xff0c;一名员工撰写一QA文档&#xff0c;并没有将它…

CDH 6.3.2升级Flink到1.17.1版本

CDH&#xff1a;6.3.2 原来的Flink&#xff1a;1.12 要升级的Flink&#xff1a;1.17.1 操作系统&#xff1a;CentOS Linux 7 一、Flink1.17编译 build.sh文件&#xff1a; #!/bin/bash set -x set -e set -vFLINK_URLsed /^FLINK_URL/!d;s/.*// flink-parcel.properties FLI…

【模板语法+数据绑定+el与data的两种写法+MVVM模型】

模板语法数据绑定el与data的两种写法MVVM模型 1 模板语法1.1 插值语法1.2 指令语法 2 数据绑定2.1 单向数据绑定2.2 双向数据绑定 3 el与data的两种写法4 MVVM模型 1 模板语法 1.1 插值语法 双大括号表达式功能&#xff1a;用于解析标签体内容语法&#xff1a;{{xxx}}&#x…

Postman的高级用法—Runner的使用

1.首先在postman新建要批量运行的接口文件夹&#xff0c;新建一个接口&#xff0c;并设置好全局变量。 2.然后在Test里面设置好要断言的方法 如&#xff1a; tests["Status code is 200"] responseCode.code 200; tests["Response time is less than 10000…

分子相互作用的人工智能

8 分子相互作用的人工智能 正如在第 5、6 和 7 节中所描述的那样&#xff0c;人工智能已经彻底改变了分子学习、蛋白质科学和材料科学领域。尽管已经广泛研究了用于单个分子的人工智能&#xff0c;但分子的物理和生物功能通常是由它们与其他分子的相互作用驱动的。在本节中&am…

iPhone数据丢失怎么办?9 佳免费 iPhone 数据恢复软件可收藏

您是否知道有多种原因可能导致 iPhone 上存储的数据永久丢失&#xff1f;然而&#xff0c;使用一些最好的免费 iPhone 数据恢复软件&#xff0c;您仍然可以恢复它。 由于我们几乎总是保存手机上的所有内容&#xff08;从联系人到媒体文件&#xff09;&#xff0c;因此 iPhone …

MySQL学习笔记27

MySQL主从复制的核心思路&#xff1a; 1、slave必须安装相同版本的mysql数据库软件。 2、master端必须开启二进制日志&#xff0c;slave端必须开启relay log 日志。 3、master主服务器和slave从服务器的server-id号不能一致。 4、slave端配置向master端来同步数据。 master…

【Java 进阶篇】MySQL 数据控制语言(DCL):管理用户权限

MySQL 是一个强大的关系型数据库管理系统&#xff0c;提供了丰富的功能和选项来管理数据库和用户。数据库管理员&#xff08;DBA&#xff09;通常使用数据控制语言&#xff08;Data Control Language&#xff0c;简称 DCL&#xff09;来管理用户的权限和访问。 本文将详细介绍…

基于php+MySql实现简易留言板(超级详细!)

基于phpMySql实现简易留言板&#xff08;超级详细&#xff01;&#xff09; 这篇文章是基于PHP和MySQL实现的一个简易留言板。该留言板具有用户注册、登录、发表留言以及查看留言等功能。首先&#xff0c;用户可以通过注册功能创建自己的账号&#xff0c;然后使用该账号进行登…

WPF 实现点击按钮跳转页面功能

方法1. 配置环境 首先添加prism依赖项&#xff0c;配置好所有文件。需要配置的有两个文件&#xff1a;App.xaml.cs和App.xaml App.xaml.cs using System.Data; using System.Linq; using System.Threading.Tasks; using System.Windows;namespace PrismDemo {/// <summa…

Linux 集锦 之 最常用的几个命令

Linux最常用的几个命令 ​ Linux系统中的命令那是相当地丰富&#xff0c;不同的版本可能还有不同的命令&#xff0c;不过Linux核心自带的命令大概有几百个&#xff0c;这个不管是什么发行版一般都是共用的。 ​ 如果希望探索Linux的所有命令&#xff0c;可能不太实际&#xf…