Hi，
I have expression (TPM) data and genotype data (in plink format) for 209 samples. I want to conduct eQTL analysis using fastlmm. The genotype files are these three: genotype_data_qc.bed, genotype_data_qc.bim, and genotype_data_qc.fam. The expression file is normalized_expression.csv, and all have been filtered. Here is my Python code:
import pandas as pd
from pysnptools.snpreader import Bed
from fastlmm.association import single_snp
import os
import numpy as np
genotype_file = Bed("genotype_data_qc.bed", count_A1=False).read()
expression_file = "normalized_expression.csv"
phenotypes = pd.read_csv(expression_file, sep="\t", index_col=0)
genotype_samples = genotype_file.iid[:, 1] # .fam 文件中的 IID
phenotype_samples = phenotypes.index
common_samples = np.intersect1d(genotype_samples, phenotype_samples)
if len(common_samples) == 0:
print("Error: No overlapping samples between genotype and phenotype.")
print("Genotype samples:", genotype_samples)
print("Phenotype samples:", phenotype_samples)
exit()
if len(common_samples) < 10:
print(f"Error: Only {len(common_samples)} samples in common, which is insufficient for analysis.")
exit()
output_dir = "eQTL_results"
os.makedirs(output_dir, exist_ok=True)
for gene in phenotypes.columns:
pheno = phenotypes[gene]
if pheno.isnull().any():
print(f"Warning: Missing values in phenotype for gene {gene}. Removing samples with missing values.")
pheno = pheno.dropna()
pheno_common = pheno.loc[common_samples]
genotype_common = genotype_file[genotype_file.iid_to_index(common_samples)]
try:
results = single_snp(
test_snps=genotype_common, # 基因型数据
pheno=pheno_common, # 表型数据
output_file_name=os.path.join(output_dir, f"eQTL_results_{gene}.txt") # 输出文件
)
print(f"eQTL analysis completed for gene: {gene}")
except Exception as e:
print(f"Error running eQTL analysis for gene {gene}: {e}")
Hi，
I have expression (TPM) data and genotype data (in plink format) for 209 samples. I want to conduct eQTL analysis using fastlmm. The genotype files are these three: genotype_data_qc.bed, genotype_data_qc.bim, and genotype_data_qc.fam. The expression file is normalized_expression.csv, and all have been filtered. Here is my Python code:
import pandas as pd
from pysnptools.snpreader import Bed
from fastlmm.association import single_snp
import os
import numpy as np
genotype_file = Bed("genotype_data_qc.bed", count_A1=False).read()
expression_file = "normalized_expression.csv"
phenotypes = pd.read_csv(expression_file, sep="\t", index_col=0)
genotype_samples = genotype_file.iid[:, 1] # .fam 文件中的 IID
phenotype_samples = phenotypes.index
common_samples = np.intersect1d(genotype_samples, phenotype_samples)
if len(common_samples) == 0:
print("Error: No overlapping samples between genotype and phenotype.")
print("Genotype samples:", genotype_samples)
print("Phenotype samples:", phenotype_samples)
exit()
if len(common_samples) < 10:
print(f"Error: Only {len(common_samples)} samples in common, which is insufficient for analysis.")
exit()
output_dir = "eQTL_results"
os.makedirs(output_dir, exist_ok=True)
for gene in phenotypes.columns:
pheno = phenotypes[gene]