You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

基于基因权重文件绘制带参考值红线的基因组织特异性直方图遇到问题

基于基因权重文件绘制带参考值红线的基因组织特异性直方图遇到问题

我在处理直方图代码时遇到了问题,抱歉问题可能比较基础,我是新手。

我有这样一个表格:

基因参考值表格

还有一个文件目录,里面的文件命名格式是这样的:en_Brain_Amygdala.db_brain_predict_1000.txt_european_brain_predict_1000(我知道这个命名很乱,是之前另一段代码的问题)。这些文件有两列:一列是基因名称,一列是权重数值。

我想要实现的是:用这些文件创建按组织/基因分开的直方图(每个组织的每个基因单独做直方图),用权重数值绘制直方图,并且在直方图上添加一条红色竖线,对应表格里的参考值。

下面是我目前写的代码:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from concurrent.futures import ProcessPoolExecutor

# Parameters and directories
Pnr = 2
part_dir = f"/mnt/hdd_3tb/GenomeData/Participant00{Pnr}/Results/1000genomes_PrediXcan/"
output_dir = f"/mnt/hdd_3tb/GenomeData/Participant00{Pnr}/Pictures/"
df_genes = pd.read_csv('/mnt/internserver1_data/ReportDatabase/UniProt/gene_list.csv', sep=',')

# Extract SCN-related genes
df_ENSG = df_genes[df_genes['gene'].str.contains('SCN', na=False)]
df_ENSG_dict = df_ENSG.set_index('gene')['ENSG'].to_dict()

# Match prediction and header files
brain_dbs = [
    f for f in os.listdir(part_dir)
    if f.startswith("en_Brain") and f.endswith("european_brain_predict_1000")
]

header_dbs = [
    f for f in os.listdir(part_dir)
    if f.startswith("en_Brain") and f.endswith("brain_predict_1000.txt")
]

# Create a mapping from prediction files to their corresponding header files
def normalize_name(filename):
    return filename.replace("_european_brain_predict_1000", "").replace(".db_brain_predict_1000.txt", "").replace(".txt", "")

# Create a mapping from header base names to full filenames
header_dict = {
    normalize_name(header): header for header in header_dbs
}


# Load the dataframe for the histograms with red lines
dataframe_path = f"/mnt/hdd_3tb/GenomeData/Participant00{Pnr}/Results/filtered_df_all.csv"
df_all = pd.read_csv(dataframe_path, index_col=0)
df_all.columns = [col.split('.')[0] for col in df_all.columns]  # Clean gene names (remove version)
# Function to generate histogram for each file
def generate_histogram(db_file,df, output_dir, df_all, df_ENSG_dict):
   db_base_name = db_file.replace("_european_brain_predict_1000", "")
   red_line_plotted_for_file = False
    
   for gene, ENSG in df_ENSG_dict.items():
        
        # Check if ENSG is part of any column name in the dataframe (as a substring match)
        matching_columns = [col for col in df.columns if ENSG[:10] in col]  # Partial matching using the first 10 characters of the ENSG ID
        
        if matching_columns:
            for col in matching_columns:
                column_data = df[col]
          
                # Create the histogram plot and save it directly if the gene matches
                fig, ax = plt.subplots()

                ax.hist(column_data)
                ax.set_xlabel('Predicted Gene Expression')
                ax.set_ylabel('Number of People')

                # Check if we need to add a red line
                if gene in df_all.columns:
                    tissue_name_from_file = db_base_name.split('_')[1]
                    matched_tissues = [tissue for tissue in df_all.index if tissue_name_from_file in tissue]
                    # Debugging: Check the tissues being matched
                    if matched_tissues:
                        matched_tissue = matched_tissues[0]
                        value = df_all.loc[matched_tissue, gene]
                        
                        try:
                            value = float(value)
                            if not pd.isna(value):
                                ax.axvline(x=value, color='red', linestyle='--', linewidth=2, label=f'Value: {value}')
                                ax.legend()
                                red_line_plotted_for_file = True
                        except ValueError:
                            print(f"Skipping non-numeric value for gene {gene} in {db_base_name}")
                            pass

                # Debug print to ensure saving
                if red_line_plotted_for_file:
                    output_path = os.path.join(output_dir, f"{db_base_name}_{col}_{gene}_histogram.png")
                    fig.savefig(output_path, bbox_inches='tight', pad_inches=0)
                    plt.close(fig)  # Close the figure after saving
                else:
                    print(f"No red line plotted for {gene} in {db_base_name}, skipping save.")

# Function to read and process each file
def process_files(db_file):
    db_base_name = normalize_name(db_file)
    header_filename = header_dict.get(db_base_name)
    
    if header_filename:
        header_file_path = os.path.join(part_dir, header_filename)
        if os.path.exists(header_file_path):
           df_header = pd.read_csv(header_file_path, sep='\t', nrows=0)
           column_names = df_header.columns
        else:
            print(f"Header file {header_file_path} does not exist!")
            return
    else:
        print(f"Header file for {db_file} not found in header_dict!")
        return

    # Read the prediction file
    data_file_path = os.path.join(part_dir, db_file)
    df = pd.read_csv(data_file_path, sep='\t', header=None)
    df.columns = column_names  # Apply headers

    # Call the function to generate histograms
    generate_histogram(db_file, df, output_dir, df_all, df_ENSG_dict)


# Use parallel processing to speed up file processing
with ProcessPoolExecutor() as executor:
   executor.map(process_files, brain_dbs)

print("Histograms with red lines have been generated and saved.")

备注:内容来源于stack exchange,提问作者Evangelia146

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.04.14 13:58:00