第一部分

1.首先第一步我獲得了所有下機數(shù)據(jù)的內源基因組覆蓋度和測序深度，利用bamdst軟件對bam文件進行計算

2.第二步我利用一個python腳本計算了所有質控后下機數(shù)據(jù)的數(shù)據(jù)量（bp）和與參考基因組相比（做除法）的測序深度，我明白了下機數(shù)據(jù)的測序深度就是將正反兩條序列的數(shù)據(jù)量相加再除以參考基因組的數(shù)據(jù)量。所寫的腳本如下，使用nohup python xxx.py & 運行后結果保存在nohup.out中。

import gzip

def count_bases(fastq_file):
    base_count = 0
    with gzip.open(fastq_file, "rt") as f:
        for i, line in enumerate(f):
            # Only process sequence lines (1-based: lines 2, 6, 10, ...)
            if i % 4 == 1:
                base_count += len(line.strip())
    return base_count

# 文件路徑
file1 = "dbh_1761_1P.fq.gz"
file2 = "dbh_1761_2P.fq.gz"

# 統(tǒng)計堿基數(shù)
bases_file1 = count_bases(file1)
bases_file2 = count_bases(file2)
total_bases = bases_file1 + bases_file2

# 輸出結果
print(f"File {file1} bases: {bases_file1}")
print(f"File {file2} bases: {bases_file2}")
print(f"Total bases: {total_bases}")

# 測序深度計算
reference_genome_size = 2475726961  # 東北虎參考基因組大?。▎挝唬篵p）
depth = total_bases / reference_genome_size
print(f"Sequencing depth: {depth:.2f}X")

當前結果展示如下：

image.png

為了控制糞便基因組測序數(shù)據(jù)與血液數(shù)據(jù)下機數(shù)據(jù)量相同，我們利用seqtk軟件進行剪切。代碼為：

seqtk sample -s100 dbh_1761_1P.fq.gz 0.6156 | gzip > dbh_1761_1P_subsampled.fq.gz
seqtk sample -s100 dbh_1761_2P.fq.gz 0.6156 | gzip > dbh_1761_2P_subsampled.fq.gz

為了完成我的首要的兩個目標，我對剪切后的數(shù)據(jù)進行了處理。常規(guī)流程，不予展示。
得到剪切后數(shù)據(jù)的bam文件計算了基因組深度覆蓋度之后，對歸一化后的數(shù)據(jù)進行比對，主要比對糞便數(shù)據(jù)的10X基因組覆蓋度與血液數(shù)據(jù)的10X基因組覆蓋度。數(shù)據(jù)保存在step1.xlsx文件中。根據(jù)數(shù)據(jù)繪制配對條形圖：

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取Excel數(shù)據(jù)
input_file = "step1.xlsx"
df = pd.read_excel(input_file)

# 數(shù)據(jù)處理，轉換為配對形式
paired_data = df.pivot(index="ID", columns="Style", values="10X_Coverage").reset_index()

# 繪圖
# 設置寬度和位置
x = np.arange(len(paired_data))  # 樣本編號索引
width = 0.35  # 條形寬度

# 創(chuàng)建條形圖
plt.figure(figsize=(8, 6))
plt.bar(x - width / 2, paired_data["blood"], width, label="Blood", color="#7e9d31")
plt.bar(x + width / 2, paired_data["faeces"], width, label="Faeces", color="#ce8862")

# 添加樣本編號為橫坐標標簽
plt.xticks(x, paired_data["ID"], rotation=45, fontsize=10)

# 設置圖形標題和標簽
plt.ylabel("10X Coverage", fontsize=12)
plt.xlabel("Sample ID", fontsize=12)
plt.title("Paired Bar Chart of 10X Coverage (Blood vs Faeces)", fontsize=14)
plt.legend(fontsize=12)

# 添加網(wǎng)格線
plt.grid(axis="y", linestyle="--", alpha=0.7)

# 調整布局
plt.tight_layout()

# 保存為PDF文件（矢量化文字）
output_file = "step1_editable.pdf"
plt.savefig(output_file, format="pdf", bbox_inches="tight")
plt.show()

print(f"圖表已成功保存為PDF文件：{output_file}")

image.png

原始數(shù)據(jù)的bam文件計算基因組覆蓋度，計算原始數(shù)據(jù)的測序深度，比對糞便數(shù)據(jù)和基因組數(shù)據(jù)的差別，同時增加了糞便熒光定量CT值

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import rcParams

# 設置字體矢量化，便于AI軟件調整文字
rcParams['pdf.fonttype'] = 42  # 確保PDF中嵌入矢量字體
rcParams['ps.fonttype'] = 42

# 讀取Excel數(shù)據(jù)
input_file = "step2.xlsx"  # 輸入文件
df = pd.read_excel(input_file)

# 數(shù)據(jù)處理
blood_data = df[df["style"] == "blood"]
faeces_data = df[df["style"] == "faeces"]

x = np.arange(len(blood_data))  # 樣本編號索引
width = 0.35  # 條形寬度

# 創(chuàng)建圖形和坐標軸
fig, ax1 = plt.subplots(figsize=(8, 6))

# 繪制條形圖 - 10X覆蓋度
ax1.bar(x - width / 2, blood_data["10X_Coverage(%)"], width, label="Blood - 10X Coverage (%)", color="#7e9d31", alpha=0.8)
ax1.bar(x + width / 2, faeces_data["10X_Coverage(%)"], width, label="Faeces - 10X Coverage (%)", color="#ce8862", alpha=0.6)

# 繪制條形圖 - 測序深度
ax1.bar(x - width / 2, blood_data["Seq_depth(X)"], width, label="Blood - Seq Depth", color="#ced38b", alpha=0.8)
ax1.bar(x + width / 2, faeces_data["Seq_depth(X)"], width, label="Faeces - Seq Depth", color="#bca788", alpha=0.6)

# 設置左側坐標軸
ax1.set_ylabel("10X Coverage (%) & Seq Depth", fontsize=12)
ax1.set_xlabel("Sample ID", fontsize=12)
ax1.set_xticks(x)
ax1.set_xticklabels(blood_data["ID"], rotation=45, fontsize=10)
ax1.set_ylim(0, 100)  # 限制覆蓋度的最大值為100
ax1.legend(loc="upper left", fontsize=10)
ax1.grid(axis="y", linestyle="--", alpha=0.7)

# 繪制次坐標軸 - CT值
ax2 = ax1.twinx()
ax2.plot(x, faeces_data["CT_value"], 'o-', color="#775a24", label="Faeces - CT Value")
ax2.set_ylabel("CT Value", fontsize=12, color="#775a24")
ax2.tick_params(axis="y", labelcolor="#775a24")
for i, ct in enumerate(faeces_data["CT_value"]):
    if not pd.isnull(ct):
        ax2.text(x[i], ct + 0.5, f"{ct:.2f}", ha="center", color="#775a24", fontsize=10)

# 設置次坐標軸圖例
ax2.legend(loc="upper right", fontsize=10)

# 添加標題
plt.title("10X Coverage, Seq Depth, and Faeces CT Value", fontsize=14)

# 調整布局并保存為PDF
output_file = "step2.pdf"
fig.tight_layout()
plt.savefig(output_file, format="pdf", bbox_inches="tight")
plt.show()

print(f"圖表已成功保存為PDF文件：{output_file}")

image.png

第二部分

對得到的bam文件繼續(xù)跑后續(xù)的步驟直到得到VCF文件主要步驟包括如下：

nohup gatk --java-options "-Xmx100g" CombineGVCFs -R ../../ref/NEW_WDS.genome.fa -V blood.list -O AT_blood.g.vcf.gz >log1.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" CombineGVCFs -R ../../ref/NEW_WDS.genome.fa -V blood.list -O AT_blood.g.vcf.gz >log1.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" GenotypeGVCFs -R ../../ref/NEW_WDS.genome.fa -V AT_blood.g.vcf.gz -O AT_blood.call.vcf.gz --tmp-dir ./tmp >log3.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" GenotypeGVCFs -R ../../ref/NEW_WDS.genome.fa -V AT_faece.g.vcf.gz -O AT_faece.call.vcf.gz --tmp-dir ./tmp >log4.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" SelectVariants -R ../../ref/NEW_WDS.genome.fa -V AT_blood.call.vcf.gz --select-type-to-include SNP -O AT_blood.call.snp.vcf.gz --tmp-dir ./tmp >log6.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" SelectVariants -R ../../ref/NEW_WDS.genome.fa -V AT_faece.call.vcf.gz --select-type-to-include SNP -O AT_faece.call.snp.vcf.gz --tmp-dir ./tmp >log7.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" VariantFiltration -R ../../ref/NEW_WDS.genome.fa -V AT_faece.call.snp.vcf.gz --filter-expression "QD<2.0 || FS>60.0 || MQ<40.0 || ReadPosRankSum<-8.0 || MQRankSum < -12.5" --filter-name "Filter" -O AT_faece.call.snp.filter.vcf.gz --tmp-dir ./tmp >log8.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" VariantFiltration -R ../../ref/NEW_WDS.genome.fa -V AT_blood.call.snp.vcf.gz --filter-expression "QD<2.0 || FS>60.0 || MQ<40.0 || ReadPosRankSum<-8.0 || MQRankSum < -12.5" --filter-name "Filter" -O AT_blood.call.snp.filter.vcf.gz --tmp-dir ./tmp >log9.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" SelectVariants -R ../../ref/NEW_WDS.genome.fa -V AT_blood.call.snp.filter.vcf.gz --exclude-filtered -O AT_blood.call.snp.filtered.vcf.gz --tmp-dir ./tmp >log10.txt 2>&1 &
nohup gatk --java-options "-Xmx100g" SelectVariants -R ../../ref/NEW_WDS.genome.fa -V AT_faece.call.snp.filter.vcf.gz --exclude-filtered -O AT_faece.call.snp.filtered.vcf.gz --tmp-dir ./tmp >log11.txt 2>&1 &

nohup bcftools view -e 'INFO/DP > 1182 || INFO/DP < 32' AT_blood.call.snp.filtered.2allelic.vcf.gz -O z -o AT_blood.call.snp.filtered.2allelic.DP.vcf.gz >log1.txt 2>&1 &
nohup bcftools view -e 'INFO/DP > 1182 || INFO/DP < 32' AT_faece.call.snp.filtered.2allelic.vcf.gz -O z -o AT_faece.call.snp.filtered.2allelic.DP.vcf.gz >log2.txt 2>&1 &

vcftools --gzvcf AT_faece.call.snp.filtered.2allelic.DP.vcf.gz --max-missing 0.98 --recode --recode-INFO-all --out AT_faece_SNP0.02
vcftools --gzvcf AT_blood.call.snp.filtered.2allelic.DP.vcf.gz --max-missing 0.98 --recode --recode-INFO-all --out AT_blood_SNP0.02

vcftools --gzvcf AT_blood_SNP0.02.recode.vcf.gz --maf 0.05 --recode --recode-INFO-all --out AT_blood_maf
vcftools --gzvcf AT_faece_SNP0.02.recode.vcf.gz --maf 0.05 --recode --recode-INFO-all --out AT_faece_maf

現(xiàn)統(tǒng)計糞便和血液數(shù)據(jù)SNP共有SNP數(shù)量（共有SNP包括POS，REF，ALT都一致的位點）、ALT不同位點（POS、REF相同，但ALT不同）、糞便獨有SNP位點、血液獨有SNP位點。

1.共有SNP位點和ALT不一致SNP位點

## 共有SNP位點
vcftools --gzvcf AT_faece_MAF.recode.vcf.gz --positions AT_blood_maf.recode.vcf.gz --out shared --recode
### 直接得到共有SNP位點數(shù)量
bcftools view -H -v shared.recode.vcf | wc -l
### 我同樣想得到ALT不一致的位點數(shù)量，首先，提取出共有SNP位點的染色體和物理坐標，
bcftools query -f'%CHROM\t%POS\n'  shared.recode.vcf> loci.txt
其次，使用shell腳本，將糞便和血液總VCF文件中的這些位點提取出來各自組成一個新的vcf文件，
最后使用python腳本進行比對分析得到數(shù)量，同時我也將SNP的染色體和位置信息提取出來

shell腳本如下,這是提取糞便數(shù)據(jù)的

#!/bin/bash

# 定義文件名
VCF_FILE="AT_faece_MAF.recode.vcf.gz"
LOCI_FILE="loci.txt"
OUTPUT_VCF="AT_feace.vcf"

# 創(chuàng)建臨時文件，用于存儲VCF文件中的坐標信息
TEMP_FILE="temp_snps.txt"

# 提取loci.txt文件中的位置信息
awk '{print $1, $2}' OFS="\t" $LOCI_FILE > $TEMP_FILE

# 使用bcftools提取VCF文件中的相關位點
bcftools view -R $TEMP_FILE -o $OUTPUT_VCF -O v $VCF_FILE

# 刪除臨時文件
rm $TEMP_FILE

echo "提取完成，結果保存到 $OUTPUT_VCF"

比對的python腳本如下

import pysam

# 打開兩個 VCF 文件
vcf_blood = pysam.VariantFile("AT_blood.vcf")
vcf_feace = pysam.VariantFile("AT_faece.vcf")

# 存儲 feace VCF 的信息
feace_records = {}
for record in vcf_feace:
    key = (record.chrom, record.pos, record.ref)
    feace_records[key] = record.alts

# 比較血液 VCF 和糞便 VCF
discordant_count = 0
discordant_positions = []  # 用于存儲差異位點的染色體和位置

for record in vcf_blood:
    key = (record.chrom, record.pos, record.ref)
    if key in feace_records and record.alts != feace_records[key]:
        discordant_count += 1
        discordant_positions.append((record.chrom, record.pos))

# 輸出差異位點數(shù)量
print(f"Number of discordant ALT positions: {discordant_count}")

# 將差異位點的染色體和位置保存到文件
discordant_file = "discordant_positions.txt"
with open(discordant_file, "w") as f:
    f.write("#CHROM\tPOS\n")  # 文件表頭
    for chrom, pos in discordant_positions:
        f.write(f"{chrom}\t{pos}\n")

print(f"Discordant positions saved to {discordant_file}")

2.血液獨有SNP位點

vcftools --gzvcf AT_blood_maf.recode.vcf.gz --exclude-positions AT_faece_MAF.recode.vcf.gz --out blood_only --recode

3.糞便獨有SNP位點

vcftools --gzvcf AT_faece_MAF.recode.vcf.gz --exclude-positions AT_blood_maf.recode.vcf.gz --out feces_only --recode

統(tǒng)計出所有數(shù)據(jù)之后進行繪圖（繪制維恩圖）

import matplotlib.pyplot as plt
from matplotlib_venn import venn2
import matplotlib as mpl

# 設置rcParams以支持字體矢量化
mpl.rcParams['pdf.fonttype'] = 42  # 矢量化字體
mpl.rcParams['ps.fonttype'] = 42  # 矢量化字體
mpl.rcParams['font.sans-serif'] = ['Arial']  # 使用常見的字體（如Arial）

# 數(shù)據(jù)
def generate_venn_diagram(output_file):
    shared_snps = 5382802
    blood_unique_snps = 214215
    feces_unique_snps = 146340

    # 創(chuàng)建維恩圖
    plt.figure(figsize=(8, 6))
    venn = venn2(subsets=(blood_unique_snps, feces_unique_snps, shared_snps),
                 set_labels=('Blood Genome', 'Feces Genome'))

    # 自定義顏色和樣式
    venn.get_patch_by_id('10').set_color('#cd8862')  # 血液基因組獨有部分
    venn.get_patch_by_id('01').set_color('#775b24')  # 糞便基因組獨有部分
    venn.get_patch_by_id('11').set_color('#7d9d33')  # 共享部分

    venn.get_patch_by_id('10').set_alpha(0.8)
    venn.get_patch_by_id('01').set_alpha(0.8)
    venn.get_patch_by_id('11').set_alpha(0.8)

    # 自定義字體和標簽位置
    venn.get_label_by_id('10').set_text(f'{blood_unique_snps}\nUnique to Blood')
    venn.get_label_by_id('01').set_text(f'{feces_unique_snps}\nUnique to Feces')
    venn.get_label_by_id('11').set_text(f'{shared_snps}\nShared SNPs')
    for label in venn.set_labels:
        label.set_fontsize(14)  # 集合名稱字體大小
    for label in venn.subset_labels:
        if label:
            label.set_fontsize(12)  # 區(qū)域內字體大小

    # 添加標題和美化
    plt.title("Comparison of SNPs between Blood and Feces Genome", fontsize=16, weight='bold')
    plt.tight_layout()

    # 導出為PDF
    plt.savefig(output_file, format='pdf')
    plt.show()
# 調用函數(shù)生成PDF文件
generate_venn_diagram("step4_venn_diagram.pdf")

print("Venn diagram saved as venn_diagram.pdf")

image.png

第三部分

我已經(jīng)對群體的數(shù)據(jù)進行了比對，現(xiàn)在需要進一步比對個體的差異，目標是比較每個個體兩組數(shù)據(jù)的靈敏度、特異性、測序錯誤率，并最終繪制小提琴圖

1.基于共有位點，根據(jù)共有位點坐標各自提取一個新的共有位點vcf文件，并進行后續(xù)拆分

## 在第二部分1中已經(jīng)介紹了如何根據(jù)位點坐標提取相應坐標的VCF文件，不再贅述

## 因為后面想要對兩個vcf文件進行對比，因此需要對vcf文件進行簡單處理，比如清除除基因型以外的信息，將基因型如0|0改為0/0
## 清除冗余信息
bcftools annotate -x ^FORMAT/GT -o AT_blood_clean.vcf -O v AT_blood.vcf
## 修改基因型（|改為/）
sed 's/|/\//g' AT_blood_clean.vcf > AT_blood_clean_fixed.vcf
## 進行個體拆分，生成個體VCF文件
bcftools view -s CXX_46 ../AT_faece_clean_fixed.vcf -o CXX_46.vcf
## 拆分完成后我使用一個python腳本，統(tǒng)計出了相同的位點，假陽性位點和假陰性位點的數(shù)量和位置
PYTHON腳本如下

import pysam

# 文件路徑
vcf1_path = "CXX_46.vcf"
vcf2_path = "AT_1367.vcf"
false_positive_output = "false_positive.txt"
false_negative_output = "false_negative.txt"

# 初始化統(tǒng)計變量
same_genotype_count = 0
false_positive_count = 0
false_negative_count = 0

# 打開輸出文件
with open(false_positive_output, "w") as fp_out, open(false_negative_output, "w") as fn_out:
    # 寫入標題行
    fp_out.write("#chrom\tpos\n")
    fn_out.write("#chrom\tpos\n")

    # 打開兩個 VCF 文件
    vcf1 = pysam.VariantFile(vcf1_path)
    vcf2 = pysam.VariantFile(vcf2_path)

    # 確保染色體和位點一致的遍歷
    for record1, record2 in zip(vcf1, vcf2):
        if record1.chrom != record2.chrom or record1.pos != record2.pos:
            raise ValueError(f"Position mismatch: {record1.chrom}:{record1.pos} != {record2.chrom}:{record2.pos}")

        # 提取基因型
        genotype1 = record1.samples[0]["GT"]
        genotype2 = record2.samples[0]["GT"]

        # 統(tǒng)計相同基因型
        if genotype1 == genotype2:
            same_genotype_count += 1

        # 檢查 false positive
        if genotype2 == (0, 0) and genotype1 in [(0, 1), (1, 1)]:
            false_positive_count += 1
            fp_out.write(f"{record1.chrom}\t{record1.pos}\n")

        # 檢查 false negative
        if genotype1 == (0, 0) and genotype2 in [(0, 1), (1, 1)]:
            false_negative_count += 1
            fn_out.write(f"{record2.chrom}\t{record2.pos}\n")

# 打印統(tǒng)計結果
print("Results:")
print(f"1. Same genotype count: {same_genotype_count}")
print(f"2. False positives (AT_1367: 0/0, CXX_46: 0/1 or 1/1): {false_positive_count}")
print(f"3. False negatives (CXX_46: 0/0, AT_1367: 0/1 or 1/1): {false_negative_count}")

2.統(tǒng)計完成后開始計算個體的三個指標

個體三個指標的計算公式

image.png

## 第一個是血液相同文件中（0/1+1/1）
import pysam

# 文件路徑
vcf_path = "AT_he24.vcf"

# 初始化計數(shù)器
genotype_0_1_count = 0
genotype_1_1_count = 0

# 打開 VCF 文件
vcf = pysam.VariantFile(vcf_path)

# 遍歷每個記錄
for record in vcf:
    # 獲取樣本的基因型（假設第一個樣本）
    genotype = record.samples[0]["GT"]
    
    # 統(tǒng)計基因型
    if genotype == (0, 1):  # 基因型為 0/1
        genotype_0_1_count += 1
    elif genotype == (1, 1):  # 基因型為 1/1
        genotype_1_1_count += 1

# 計算總數(shù)
total_count = genotype_0_1_count + genotype_1_1_count

# 輸出統(tǒng)計結果
print(f"Genotype 0/1 count: {genotype_0_1_count}")
print(f"Genotype 1/1 count: {genotype_1_1_count}")
print(f"Total 0/1 and 1/1 count: {total_count}")

### 第二個是血液單獨文件中（0/1+1/1）
import pysam

# 文件路徑
vcf_path = "AT_he24_blood.vcf"

# 初始化計數(shù)器
genotype_0_1_count = 0
genotype_1_1_count = 0

# 打開 VCF 文件
vcf = pysam.VariantFile(vcf_path)

# 遍歷每個記錄
for record in vcf:
    # 獲取樣本的基因型（假設第一個樣本）
    genotype = record.samples[0]["GT"]
    
    # 統(tǒng)計基因型
    if genotype == (0, 1):  # 基因型為 0/1
        genotype_0_1_count += 1
    elif genotype == (1, 1):  # 基因型為 1/1
        genotype_1_1_count += 1

# 計算總數(shù)
total_count = genotype_0_1_count + genotype_1_1_count

# 輸出統(tǒng)計結果
print(f"Genotype 0/1 count: {genotype_0_1_count}")
print(f"Genotype 1/1 count: {genotype_1_1_count}")
print(f"Total 0/1 and 1/1 count: {total_count}")

### 第三個是糞便單獨文件中（0/1+1/1）
import pysam

# 文件路徑
vcf_path = "H-11-1386-HL.feces.vcf"

# 初始化計數(shù)器
genotype_0_1_count = 0
genotype_1_1_count = 0

# 打開 VCF 文件
vcf = pysam.VariantFile(vcf_path)

# 遍歷每個記錄
for record in vcf:
    # 獲取樣本的基因型（假設第一個樣本）
    genotype = record.samples[0]["GT"]
    
    # 統(tǒng)計基因型
    if genotype == (0, 1):  # 基因型為 0/1
        genotype_0_1_count += 1
    elif genotype == (1, 1):  # 基因型為 1/1
        genotype_1_1_count += 1

# 計算總數(shù)
total_count = genotype_0_1_count + genotype_1_1_count

# 輸出統(tǒng)計結果
print(f"Genotype 0/1 count: {genotype_0_1_count}")
print(f"Genotype 1/1 count: {genotype_1_1_count}")
print(f"Total 0/1 and 1/1 count: {total_count}")

得到了靈敏度，差異性與測序錯誤率的數(shù)值，進行繪圖，對靈敏度和差異性我繪制了小提琴圖，對測序錯誤率我結合了10X基因組覆蓋度差值繪制了散點圖
靈敏度和差異性小提琴圖

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 確保PDF中的字體可編輯
rcParams['ps.fonttype'] = 42

# 讀取數(shù)據(jù)
input_file = "step5_3_1.xlsx"
data = pd.read_excel(input_file)

# 數(shù)據(jù)格式轉換
df_melted = data.melt(id_vars=["ID"], value_vars=["Sensitivity", "Idiosyncrasy"],
                      var_name="Metric", value_name="Value")

# 自定義配色
custom_palette = ["#7e9d31", "#775a24"]  # 藍色和橙色

# 創(chuàng)建小提琴圖
plt.figure(figsize=(8, 6))
sns.violinplot(x="Metric", y="Value", data=df_melted, inner="quartile", palette=custom_palette)

# 設置圖表標題和軸標簽
plt.title("Distribution of Sensitivity and Idiosyncrasy", fontsize=14)
plt.ylabel("Value", fontsize=12)
plt.xlabel("Metric", fontsize=12)
plt.ylim(0.94, 1.00)

# 添加基準線
plt.axhline(0.95, color="red", linestyle="--", linewidth=1, label="Threshold (0.95)")
plt.legend(fontsize=10)

# 添加網(wǎng)格
plt.grid(axis="y", linestyle="--", alpha=0.6)

# 調整布局
plt.tight_layout()

# 保存圖形為矢量化的PDF格式
output_file = "step5_3_1.pdf"
plt.savefig(output_file, format="pdf")
print(f"圖形已成功保存為 {output_file}，可在AI軟件中編輯字體和文字大小。")
plt.show()

image.png

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 讀取文件
data_file = "step5_3_2.xlsx"
data = pd.read_excel(data_file)

# 設置字體和矢量化輸出
rcParams['pdf.fonttype'] = 42  # 矢量化輸出
rcParams['ps.fonttype'] = 42
plt.rcParams['font.sans-serif'] = ['Arial']  # 設置字體為Arial

# 數(shù)據(jù)準備
highlight_ids = {"AT_1367", "AT_1578", "AT_1716", "AT_he24"}

x = data['10X genome coverage difference']
y = data['Sequencing Error Rate']
labels = data['ID']

# 配色設置
highlight_color = "#587d03"  # 高亮顏色
normal_color = "#ced38c"  # 普通點顏色

# 繪制散點圖
plt.figure(figsize=(8, 6))
for xi, yi, label in zip(x, y, labels):
    if label in highlight_ids:
        plt.scatter(xi, yi, color=highlight_color, s=100, label=label if label not in plt.gca().get_legend_handles_labels()[1] else None)
        plt.text(xi, yi, label, fontsize=9, ha='right')
    else:
        plt.scatter(xi, yi, color=normal_color, s=50)

# 添加輔助線
plt.axhline(y=y.mean(), color='gray', linestyle='--', label="Avg Error Rate")
plt.axvline(x=x.mean(), color='gray', linestyle='--', label="Avg Coverage Diff")

# 圖表設置
plt.title("Sequencing Error Rate vs 10X Genome Coverage Difference", fontsize=14)
plt.xlabel("10X Genome Coverage Difference", fontsize=12)
plt.ylabel("Sequencing Error Rate", fontsize=12)
plt.legend(title="Highlighted IDs", fontsize=10)
plt.grid(alpha=0.3)

# 保存為PDF
output_file = "step5_3_2.pdf"
plt.tight_layout()
plt.savefig(output_file, format='pdf', bbox_inches='tight')

# 顯示圖表
plt.show()

print(f"圖表已保存為 {output_file}")

image.png

第四部分

我們已經(jīng)完成了之前的深度比對以及SNP位點錯誤率等的比對，下面開始進行進一步的遺傳參數(shù)比對

1.等位基因頻率分布比對，用疊加密度曲線圖展示糞便數(shù)據(jù)和血液數(shù)據(jù)之間的異同
等位基因頻率分布是指不同頻率區(qū)間的等位基因數(shù)量分布情況，即樣本中具有不同頻率的變異等位基因在群體中的比例，該分布反映了遺傳多樣性在不同等位基因頻率下的分布特征。
其實就是統(tǒng)計次級等位基因頻率
使用vcftools軟件計算等位基因頻率

### 計算等位基因頻率
vcftools --gzvcf ../../var/AT_faece_MAF.recode.vcf.gz --freq --out faece_allele_freq

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 文件路徑
faece_file_path = "faece_allele_freq.frq"  # 糞便數(shù)據(jù)文件路徑
blood_file_path = "blood_allele_freq.frq"  # 血液數(shù)據(jù)文件路徑

# 確保輸出為矢量化字體，方便在 AI 軟件中修改字體
rcParams['pdf.fonttype'] = 42  # 嵌入矢量字體
rcParams['ps.fonttype'] = 42  # 嵌入矢量字體

# 函數(shù)：提取次級等位基因頻率 (< 0.5)
def extract_minor_allele_freq(file_path):
    df = pd.read_csv(file_path, sep="\t")  # 讀取 .frq 文件
    frequencies = []  # 存儲次級等位基因頻率
    for freq_str in df['{ALLELE:FREQ}']:
        alleles = freq_str.split()  # 切分字符串，如 "A:0.8 T:0.2"
        freqs = [float(a.split(":")[1]) for a in alleles]  # 提取頻率值
        minor_freq = min(freqs)  # 取較小頻率值，次級等位基因頻率
        if minor_freq < 0.5:  # 保留次級等位基因頻率 < 0.5 的數(shù)值
            frequencies.append(minor_freq)
    return frequencies

# 提取糞便和血液樣本的次級等位基因頻率
faece_frequencies = extract_minor_allele_freq(faece_file_path)
blood_frequencies = extract_minor_allele_freq(blood_file_path)

# 創(chuàng)建用于繪圖的數(shù)據(jù)框
combined_df = pd.DataFrame({
    'Frequency': faece_frequencies + blood_frequencies,
    'Sample Type': ['Faeces'] * len(faece_frequencies) + ['Blood'] * len(blood_frequencies)
})

# 繪制疊加密度曲線圖
plt.figure(figsize=(4, 4))

# 設置自定義配色
custom_palette = {'Faeces': '#7e9d31', 'Blood': '#ce8862'}  # 可自行調整色號

sns.kdeplot(data=combined_df, x='Frequency', hue='Sample Type', fill=True, alpha=0.6, bw_adjust=1.2, palette=custom_palette)
plt.xlabel('Minor Allele Frequency', fontsize=14)
plt.ylabel('Density', fontsize=14)
plt.title('Density Plot of Minor Allele Frequency (Faeces vs Blood)', fontsize=16)

# 去掉圖例外邊框，優(yōu)化間距
plt.legend(title='Sample Type', frameon=False, loc='upper right')
plt.tight_layout()  # 自動調整間距

# 保存圖表為 PDF 格式
plt.savefig("allele_frequency_density_plot_customized.pdf", format='pdf', bbox_inches='tight')
plt.show()

image.png

## 第一步從VCF文件中提取數(shù)據(jù)
import pysam
import pandas as pd

# 文件路徑
faece_vcf_path = "AT_faece_MAF.recode.vcf.gz"  # 糞便數(shù)據(jù)
blood_vcf_path = "AT_blood_maf.recode.vcf.gz"  # 血液數(shù)據(jù)

# 函數(shù)：計算 SNP 密度
def calculate_snp_density(vcf_path, window_size=50000):
    vcf = pysam.VariantFile(vcf_path)  # 讀取 VCF 文件
    density_data = []  # 存儲 SNP 密度數(shù)據(jù)

    for chrom in vcf.header.contigs:  # 遍歷每個染色體
        chrom_length = vcf.header.contigs[chrom].length  # 獲取染色體長度
        for start in range(0, chrom_length, window_size):
            end = min(start + window_size, chrom_length)
            snp_count = 0
            for record in vcf.fetch(chrom, start, end):  # 獲取該窗口內的 SNP 記錄
                snp_count += 1
            density_data.append([chrom, start, end, snp_count])
    
    vcf.close()
    return pd.DataFrame(density_data, columns=['Chromosome', 'Start', 'End', 'SNP_Count'])

# 計算 SNP 密度
faece_density = calculate_snp_density(faece_vcf_path, window_size=50000)
blood_density = calculate_snp_density(blood_vcf_path, window_size=50000)

# 添加樣本類型
faece_density['Sample_Type'] = 'Feces'
blood_density['Sample_Type'] = 'Blood'

# 合并數(shù)據(jù)并保存
combined_density = pd.concat([faece_density, blood_density], ignore_index=True)
combined_density.to_csv("snp_density_plot_input.csv", index=False)

print("SNP 密度數(shù)據(jù)已保存為 snp_density_plot_input.csv")

第二步，將生成文件中無用的短染色體去除

import pandas as pd

# 讀取 CSV 文件
file_path = "snp_density_plot_input.csv"
df = pd.read_csv(file_path)

# 提取符合條件的行（染色體編號為 Chrwds1, Chrwds2, Chrwds3, Chrwds4）
chromosomes_to_keep = ['Chrwds1', 'Chrwds2', 'Chrwds3', 'Chrwds4', 'Chrwds5', 'Chrwds6', 'Chrwds7', 'Chrwds8', 'Chrwds9', 'Chrwds10', 'Chrwds11', 'Chrwds12', 'Chrwds13', 'Chrwds14', 'Chrwds15', 'Chrwds16', 'Chrwds17', 'Chrwds18', 'Chrwds19', 'Chrwds20']
filtered_df = df[df['Chromosome'].isin(chromosomes_to_keep)]

# 保存提取的內容為新 CSV 文件
output_file_path = "filtered_snp_density.csv"
filtered_df.to_csv(output_file_path, index=False)

print(f"提取完成，已保存為 {output_file_path}")

第三步，繪圖

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取 CSV 文件
df = pd.read_csv("step6_2_filtered_snp_density.csv")

# 計算 SNP 密度
window_size_kb = 100  # 100 kb
df['SNP_Density'] = df['SNP_Count'] / window_size_kb  # 每 kb 的 SNP 密度

# 為每個染色體編號一個連續(xù)的起始位置索引，使 X 軸展示染色體編號
chromosome_order = [f'Chrwds{i}' for i in range(1, 21)]  # Chrwds1 到 Chrwds20
df['Chromosome'] = pd.Categorical(df['Chromosome'], categories=chromosome_order, ordered=True)

# 添加新列：每個染色體分段的位置索引
chromosome_offset = {}
offset = 0

for chrom in chromosome_order:
    chrom_length = df[df['Chromosome'] == chrom]['Start'].max()  # 每條染色體的最大起始位置
    chromosome_offset[chrom] = offset
    offset += chrom_length + 1e6  # 間隔 1 Mb

df['Relative_Position'] = df.apply(lambda row: row['Start'] + chromosome_offset[row['Chromosome']], axis=1)

# 數(shù)據(jù)平滑：滾動平均
df['SNP_Density_Smoothed'] = df.groupby('Chromosome')['SNP_Density'].transform(lambda x: x.rolling(10, min_periods=1).mean())

# 自定義配色（用色號指定）
custom_palette = {
    'Faeces': '#0f9dd0',  # 深棕色
    'Blood': '#f6c81f'    # 綠色
}

# 將糞便數(shù)據(jù)放在最頂層：調整繪圖順序
df['Sample_Type'] = pd.Categorical(df['Sample_Type'], categories=['Blood', 'Faeces'], ordered=True)

# 繪圖：用 `Relative_Position` 代替 `Start`
plt.figure(figsize=(10, 6))
sns.lineplot(
    data=df, 
    x='Relative_Position', 
    y='SNP_Density_Smoothed', 
    hue='Sample_Type', 
    ci=None, 
    linewidth=0.8, 
    palette=custom_palette, 
    alpha=0.6  # 設置透明度
)

# 設置染色體標簽為 X 軸刻度
chrom_labels = {v: k for k, v in chromosome_offset.items()}
plt.xticks(list(chrom_labels.keys()), list(chrom_labels.values()), rotation=45)
plt.xlabel('Chromosome', fontsize=14)
plt.ylabel('SNP Density (per kb)', fontsize=14)
plt.title('SNP Density Distribution by Chromosome (Faeces vs Blood)', fontsize=16)
plt.legend(title='Sample Type', loc='upper right')
plt.tight_layout()

# 保存為矢量化 PDF 格式
plt.savefig("snp_density_plot_with_alpha_overlay.pdf", format='pdf', bbox_inches='tight')
plt.show()

image.png

import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import linregress
import numpy as np
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取Excel文件數(shù)據(jù)
input_file = "step7_1HET.xlsx"  # 輸入文件名
sheet_name = 0  # 默認第一個工作表
df = pd.read_excel(input_file, sheet_name=sheet_name)

# 檢查數(shù)據(jù)列名
if not {"Blood", "Feces"}.issubset(df.columns):
    raise ValueError("Excel 文件缺少 'Blood' 和 'Feces' 列，請檢查文件格式。")

# 進行線性回歸
slope, intercept, r_value, p_value, std_err = linregress(df["Blood"], df["Feces"])

# 計算擬合曲線
x = np.linspace(df["Blood"].min(), df["Blood"].max(), 100)
y_fit = slope * x + intercept

# 自定義配色選項
scatter_color = "#7e9d31"  # 默認散點顏色（可修改為其他色號）
line_color = "#775a24"     # 默認擬合線顏色（可修改為其他色號）

# 繪圖
plt.figure(figsize=(6, 6))
plt.scatter(df["Blood"], df["Feces"], color=scatter_color)
plt.plot(x, y_fit, color=line_color)
plt.xlabel("Blood Genomic Heterozygosity")
plt.ylabel("Feces Genomic Heterozygosity")
plt.title("Linear Regression Analysis")
plt.legend()

# 設置坐標軸框線可見性
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(True)
ax.spines['bottom'].set_visible(True)

# 確?？潭染€在正確位置顯示
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

# 設置外側刻度線
ax.tick_params(axis='both', which='both', direction='out', length=6, width=1, colors='black')

# 設置科學計數(shù)法格式
plt.ticklabel_format(style='sci', axis='both', scilimits=(0, 0))

# 在圖中顯示 R2 和 p-value
plt.text(df["Blood"].min(), df["Feces"].max() * 0.99, f"R2 = {r_value**2:.4f}\nP-value = {p_value:.4e}", fontsize=12, bbox=dict(facecolor='white', alpha=0.5))

# 保存為矢量化 PDF
output_file = "step7_1HET.pdf"
plt.savefig(output_file, format='pdf', bbox_inches='tight')

plt.show()

print(f"圖表已成功保存為 {output_file}，并以矢量格式導出為 PDF。")

image.png

import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import linregress
import numpy as np
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取Excel文件數(shù)據(jù)
input_file = "step7_2Relationship.xlsx"  # 輸入文件名
sheet_name = 0  # 默認第一個工作表
df = pd.read_excel(input_file, sheet_name=sheet_name)

# 檢查數(shù)據(jù)列名
if not {"Blood", "Feces"}.issubset(df.columns):
    raise ValueError("Excel 文件缺少 'Blood' 和 'Feces' 列，請檢查文件格式。")

# 進行線性回歸
slope, intercept, r_value, p_value, std_err = linregress(df["Blood"], df["Feces"])

# 計算擬合曲線
x = np.linspace(df["Blood"].min(), df["Blood"].max(), 100)
y_fit = slope * x + intercept

# 自定義配色選項
scatter_color = "#7e9d31"  # 默認散點顏色（可修改為其他色號）
line_color = "#775a24"     # 默認擬合線顏色（可修改為其他色號）

# 繪圖
plt.figure(figsize=(6, 6))
plt.scatter(df["Blood"], df["Feces"], color=scatter_color)
plt.plot(x, y_fit, color=line_color)
plt.xlabel("Blood Genomic Kinship")
plt.ylabel("Feces Genomic Kinship")
plt.title("Linear Regression Analysis")
plt.legend()

# 設置坐標軸框線可見性
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(True)
ax.spines['bottom'].set_visible(True)

# 確保刻度線在正確位置顯示
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

# 設置外側刻度線
ax.tick_params(axis='both', which='both', direction='out', length=6, width=1, colors='black')

# 設置科學計數(shù)法格式
plt.ticklabel_format(style='sci', axis='both', scilimits=(0, 0))

# 在圖中顯示 R2 和 p-value
plt.text(df["Blood"].min(), df["Feces"].max() * 0.99, f"R2 = {r_value**2:.4f}\nP-value = {p_value:.4e}", fontsize=12, bbox=dict(facecolor='white', alpha=0.5))

# 保存為矢量化 PDF
output_file = "step7_2Relationship.pdf"
plt.savefig(output_file, format='pdf', bbox_inches='tight')

plt.show()

print(f"圖表已成功保存為 {output_file}，并以矢量格式導出為 PDF。")

image.png

import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import linregress
import numpy as np
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取Excel文件數(shù)據(jù)
input_file = "step7_3Froh.xlsx"  # 輸入文件名
sheet_name = 0  # 默認第一個工作表
df = pd.read_excel(input_file, sheet_name=sheet_name)

# 檢查數(shù)據(jù)列名
if not {"Blood", "Feces"}.issubset(df.columns):
    raise ValueError("Excel 文件缺少 'Blood' 和 'Feces' 列，請檢查文件格式。")

# 進行線性回歸
slope, intercept, r_value, p_value, std_err = linregress(df["Blood"], df["Feces"])

# 計算擬合曲線
x = np.linspace(df["Blood"].min(), df["Blood"].max(), 100)
y_fit = slope * x + intercept

# 自定義配色選項
scatter_color = "#7e9d31"  # 默認散點顏色（可修改為其他色號）
line_color = "#775a24"     # 默認擬合線顏色（可修改為其他色號）

# 繪圖
plt.figure(figsize=(6, 6))
plt.scatter(df["Blood"], df["Feces"], color=scatter_color)
plt.plot(x, y_fit, color=line_color)
plt.xlabel("Blood Genomic Froh")
plt.ylabel("Feces Genomic Froh")
plt.title("Linear Regression Analysis")
plt.legend()

# 設置坐標軸框線可見性
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(True)
ax.spines['bottom'].set_visible(True)

# 確?？潭染€在正確位置顯示
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

# 設置外側刻度線
ax.tick_params(axis='both', which='both', direction='out', length=6, width=1, colors='black')

# 設置科學計數(shù)法格式
plt.ticklabel_format(style='sci', axis='both', scilimits=(0, 0))

# 在圖中顯示 R2 和 p-value
plt.text(df["Blood"].min(), df["Feces"].max() * 0.99, f"R2 = {r_value**2:.4f}\nP-value = {p_value:.4e}", fontsize=12, bbox=dict(facecolor='white', alpha=0.5))

# 保存為矢量化 PDF
output_file = "step7_3Froh.pdf"
plt.savefig(output_file, format='pdf', bbox_inches='tight')

plt.show()

print(f"圖表已成功保存為 {output_file}，并以矢量格式導出為 PDF。")

image.png

import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import linregress
import numpy as np
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取Excel文件數(shù)據(jù)
input_file = "step7_4Load.xlsx"  # 輸入文件名
sheet_name = 0  # 默認第一個工作表
df = pd.read_excel(input_file, sheet_name=sheet_name)

# 檢查數(shù)據(jù)列名
if not {"Blood", "Feces"}.issubset(df.columns):
    raise ValueError("Excel 文件缺少 'Blood' 和 'Feces' 列，請檢查文件格式。")

# 進行線性回歸
slope, intercept, r_value, p_value, std_err = linregress(df["Blood"], df["Feces"])

# 計算擬合曲線
x = np.linspace(df["Blood"].min(), df["Blood"].max(), 100)
y_fit = slope * x + intercept

# 自定義配色選項
scatter_color = "#7e9d31"  # 默認散點顏色（可修改為其他色號）
line_color = "#775a24"     # 默認擬合線顏色（可修改為其他色號）

# 繪圖
plt.figure(figsize=(6, 6))
plt.scatter(df["Blood"], df["Feces"], color=scatter_color)
plt.plot(x, y_fit, color=line_color)
plt.xlabel("Blood Genomic Load")
plt.ylabel("Feces Genomic Load")
plt.title("Linear Regression Analysis")
plt.legend()

# 設置坐標軸框線可見性
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(True)
ax.spines['bottom'].set_visible(True)

# 確?？潭染€在正確位置顯示
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')

# 設置外側刻度線
ax.tick_params(axis='both', which='both', direction='out', length=6, width=1, colors='black')

# 設置科學計數(shù)法格式
plt.ticklabel_format(style='sci', axis='both', scilimits=(0, 0))

# 在圖中顯示 R2 和 p-value
plt.text(df["Blood"].min(), df["Feces"].max() * 0.99, f"R2 = {r_value**2:.4f}\nP-value = {p_value:.4e}", fontsize=12, bbox=dict(facecolor='white', alpha=0.5))

# 保存為矢量化 PDF
output_file = "step7_4Load.pdf"
plt.savefig(output_file, format='pdf', bbox_inches='tight')

plt.show()

print(f"圖表已成功保存為 {output_file}，并以矢量格式導出為 PDF。")

image.png

第5部分

開始進行功能基因分析，查看這些差異位點是否具有位置偏好，以及是否會影響免疫相關基因和其他關鍵基因。關注與健康或適應性相關的位點，分析可能的系統(tǒng)偏差。

第一步：先找到每個個體的糞便差異SNP位點在基因組的什么區(qū)域。比如外顯子，內含子以及調控區(qū)

##先抽出血液、糞便獨有文件中的0/1和1/1（以血液為例）
import pysam

# 文件路徑
vcf_path = "AT_1367_blood.vcf"
output_path = "blood.txt"

# 初始化計數(shù)器
genotype_0_1_count = 0
genotype_1_1_count = 0

# 打開輸出文件
with open(output_path, "w") as output_file:
    # 寫入表頭
    output_file.write("Chromosome\tPosition\tGenotype\n")
    
    # 打開 VCF 文件
    vcf = pysam.VariantFile(vcf_path)
    
    # 遍歷每個記錄
    for record in vcf:
        # 獲取樣本的基因型（假設第一個樣本）
        genotype = record.samples[0]["GT"]
        
        # 統(tǒng)計基因型并保存位點信息
        if genotype == (0, 1):  # 基因型為 0/1
            genotype_0_1_count += 1
            output_file.write(f"{record.chrom}\t{record.pos}\t0/1\n")
        elif genotype == (1, 1):  # 基因型為 1/1
            genotype_1_1_count += 1
            output_file.write(f"{record.chrom}\t{record.pos}\t1/1\n")

# 計算總數(shù)
total_count = genotype_0_1_count + genotype_1_1_count

# 輸出統(tǒng)計結果
print(f"Genotype 0/1 count: {genotype_0_1_count}")
print(f"Genotype 1/1 count: {genotype_1_1_count}")
print(f"Total 0/1 and 1/1 count: {total_count}")
print(f"Chromosome positions saved to {output_path}")

過濾掉得到的blood.txt文件中的無用信息

vi blood.txt (去掉首行）
awk '{print $1,$2}' blood.txt > blood_snp.txt
### 之前得到的共有SNP位點的vcf文件，也要去掉第一行
cat blood_snp.txt feces_snp.txt false_negative.txt false_positive.txt > position.txt

轉換為bed文件，進行比對

## position.txt 文件轉換為bed文件
awk '{print $1"\t"$2-1"\t"$2}' position.txt > positions.bed
## gff文件轉換為bed文件
awk '{ if ($3 == "mRNA" || $3 == "CDS") print $1"\t"$4-1"\t"$5"\t"$3"\t"$9 }' ../../../../ref/NEW_WDS.gff > NEW_WDS.bed
### 調控區(qū)域轉換為bed文件
awk '{ if ($3 == "mRNA") print $1, $4-2000, $4, "promoter_region" }' ../../../../ref/NEW_WDS.gff > promoter_regions.bed
## 將列與列之間的空格改為制表符TAB
sed -i 's/ \+/\t/g' promoter_regions.bed
## 首先比對CDS外顯子區(qū)和內含子區(qū)
bedtools intersect -a positions.bed -b NEW_WDS.bed -wa -wb > result.bed
## 統(tǒng)計數(shù)量
grep "CDS" result.bed | wc -l
## 其次比對調控區(qū)
bedtools intersect -a positions.bed -b promoter_regions.bed -wa -wb > result1.bed

統(tǒng)計完成這些每個個體的CDS數(shù)量、內含子數(shù)量以及調控區(qū)數(shù)量后，繪制了一個環(huán)形圖

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取 Excel 文件
data_file = "step8_1exon.xlsx"
df = pd.read_excel(data_file)

# 檢查數(shù)據(jù)格式是否正確
#print("Data preview:")
#print(df.head())

# 計算總和
df_sum = pd.DataFrame({
    "region": ["exon", "intron", "regulatory_region"],
    "SNP_count": [df["exon"].sum(), df["introns"].sum(), df["regulatory region"].sum()]
})

# 繪制環(huán)形圖函數(shù)
def draw_donut_chart(dataframe, colors=None, output_file="step8_1_SNP_distribution.pdf"):
    plt.figure(figsize=(8, 8))
    
    # 畫環(huán)形圖，矢量輸出
    plt.pie(
        dataframe["SNP_count"],
        labels=dataframe["region"],
        autopct='%1.1f%%',
        startangle=90,
        wedgeprops={'width': 0.4, 'edgecolor': 'w'},
        colors=colors  # 顏色可以由用戶輸入
    )
    plt.title("Distribution of SNP Counts in Different Genomic Regions", pad=20)
    
    # 保存為矢量化 PDF
    plt.savefig(output_file, format='pdf', bbox_inches='tight')
    plt.close()
    print(f"Chart saved as {output_file}")

# 示例顏色輸入
custom_colors = ["#fb876e", "#b4cd58", "#7dd8cf"]  # 用戶可以根據(jù)需要修改色號

# 調用繪制函數(shù)
draw_donut_chart(df_sum, colors=custom_colors)
plt.show()

image.png

grep "CDS" result.bed > CDS.bed
### 提取基因名
awk '{print $8}' CDS.bed | awk -F "=" '{print $2}' > gene_id.txt
## 去除重復
sort gene_id.txt | uniq > AT_1716_gene_id_move.txt
### 進行匯總后，生成一個excel表格，然后根據(jù)表格繪制了個體間相同基因數(shù)以及不同樣本頻率下，出現(xiàn)的相同基因的次數(shù)

首先，我利用統(tǒng)計的step9_1gene.xlsx表格數(shù)據(jù)生成了一個個體間相同基因矩陣csv文件，腳本如下

import pandas as pd

# 讀取 Excel 文件
file_path = "step9_1gene.xlsx"
df = pd.read_excel(file_path)

# 提取每列的基因集合，包括第一列
gene_sets = {}
for col in df.columns:  # 從第一列開始
    gene_sets[col] = set(df[col].dropna().tolist())  # 去除空值并轉為集合

# 構建交集矩陣
sample_names = list(gene_sets.keys())
intersection_matrix = pd.DataFrame(index=sample_names, columns=sample_names)

# 統(tǒng)計樣本間共享基因數(shù)量
for s1 in sample_names:
    for s2 in sample_names:
        intersection_matrix.loc[s1, s2] = len(gene_sets[s1] & gene_sets[s2])  # 取交集長度

# 保存矩陣為 Excel 文件
output_file = "gene_intersection_matrix.xlsx"
intersection_matrix.to_excel(output_file)

# 輸出矩陣文件路徑
print(f"Intersection matrix saved to: {output_file}")

接著，我利用矩陣繪制了一個環(huán)形熱圖

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import Normalize
from math import pi
from matplotlib import rcParams

# 設置字體矢量化，便于AI軟件調整文字
rcParams['pdf.fonttype'] = 42  # 確保PDF中嵌入矢量字體
rcParams['ps.fonttype'] = 42
# 讀取 Excel 文件
file_name = "step9_2_cir_gene.xlsx"
df = pd.read_excel(file_name, index_col=0)

# 獲取樣本標簽和數(shù)據(jù)矩陣
labels = df.columns
size = len(labels)
values = df.values

# 設置顏色映射范圍（可自定義）
def plot_open_circular_heatmap(df, cmap="Blues", output_file="step9_1_circular_heatmap.pdf"):
    """
    繪制不閉合環(huán)形熱圖并導出為PDF
    Args:
        df (pd.DataFrame): 輸入數(shù)據(jù)矩陣
        cmap (str): 顏色映射方案，可自定義
        output_file (str): 導出文件名
    """
    # 歸一化顏色映射
    norm = Normalize(vmin=np.min(df.values), vmax=np.max(df.values))
    colormap = plt.get_cmap(cmap)

    # 創(chuàng)建圖形
    fig, ax = plt.subplots(figsize=(12, 12), subplot_kw={'projection': 'polar'})

    # 設置環(huán)形范圍：開口角度
    opening_angle = pi / 4  # 45 度開口
    theta_start = opening_angle / 2  # 起始角度
    theta_end = 2 * pi - theta_start  # 結束角度

    # 繪制環(huán)形熱圖
    for i in range(size):
        for j in range(size):
            theta_1 = theta_start + (theta_end - theta_start) * i / size
            theta_2 = theta_start + (theta_end - theta_start) * (i + 1) / size
            r_start = 0.8 + j / size  # 留出較大的中心空白
            r_end = r_start + 1 / size

            color = colormap(norm(df.values[i, j]))
            ax.barh(r_start, theta_2 - theta_1, left=theta_1, height=r_end - r_start, color=color, edgecolor='white')

    # 設置圖形參數(shù)
    ax.set_xticks(np.linspace(theta_start, theta_end, size))
    ax.set_xticklabels(labels, fontsize=10)
    ax.set_yticks([])  # 去除虛線刻度
    ax.spines['polar'].set_visible(False)  # 去掉極坐標外框線
    ax.set_title("Open Circular Heatmap of Gene Intersection", pad=20)

    # 添加顏色條
    sm = plt.cm.ScalarMappable(cmap=colormap, norm=norm)
    cbar = fig.colorbar(sm, ax=ax, shrink=0.8, pad=0.1)
    cbar.set_label('Shared Gene Count')

    # 保存 PDF 矢量格式
    plt.savefig(output_file, format='pdf', bbox_inches='tight')
    plt.show()

# 執(zhí)行繪圖函數(shù)
plot_open_circular_heatmap(df, cmap="Blues", output_file="step9_1_circular_heatmap_gene_intersection.pdf")

接著我又統(tǒng)計了不同樣本頻次的相同基因柱狀圖

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 設置字體矢量化，便于AI軟件調整文字
rcParams['pdf.fonttype'] = 42  # 確保PDF中嵌入矢量字體
rcParams['ps.fonttype'] = 42
# 讀取 Excel 文件
file_path = "step9_1gene.xlsx"
df = pd.read_excel(file_path)

# 提取每列的基因集合，包括第一列
gene_sets = {}
for col in df.columns:  # 從第一列開始
    gene_sets[col] = set(df[col].dropna().tolist())  # 去除空值并轉為集合

# 統(tǒng)計基因在不同樣本數(shù)中出現(xiàn)的頻次
gene_presence_count = {}
for sample_set in gene_sets.values():
    for gene in sample_set:
        gene_presence_count[gene] = gene_presence_count.get(gene, 0) + 1

# 統(tǒng)計每個頻次的基因數(shù)量
frequency_distribution = pd.Series(list(gene_presence_count.values())).value_counts().sort_index()

# 將結果保存為 CSV 文件
output_file = "gene_frequency_distribution.csv"
frequency_distribution.to_csv(output_file, header=["Number of Genes"], index_label="Number of Samples")

# 繪制柱狀圖，設置柱子緊密排列
plt.figure(figsize=(6, 4))
plt.bar(frequency_distribution.index, frequency_distribution.values, color="cornflowerblue", width=0.9)  # 寬度設置為 1.0
plt.xlabel("Number of Samples Containing the Gene")
plt.ylabel("Number of Genes")
plt.title("Gene Frequency Distribution Across Samples")
plt.xticks(range(1, len(gene_sets) + 1))
barplot_file = "step9_2gene_frequency_barplot_tight.pdf"
plt.savefig(barplot_file, format='pdf', bbox_inches='tight')
plt.show()

# 輸出結果文件路徑
print(f"Frequency distribution saved to: {output_file}")
print(f"Bar plot saved to: {barplot_file}")

最后我將圖拼接到一起，生成了一張組合圖

image.png

cat AT_1367_gene_id_move.txt AT_1578_gene_id_move.txt AT_1716_gene_id_move.txt AT_1731_gene_id_move.txt AT_hd26_gene_id_move.txt AT_hd32_gene_id_move.txt AT_hd36_gene_id_move.txt AT_hd70_gene_id_move.txt AT_hd88_gene_id_move.txt AT_he24_gene_id_move.txt dbh_1761_gene_id_move.txt > gene_id.txt

## 去除重復
sort gene_id.txt | uniq > gene_id_move.txt

### 我有一個文件是基因組的完整的注釋文件，里面包含go分析和KEGG分析的基因功能和通路，
##然后我可以將上一步生成文件中的所有基因從注釋文件中提取出來，
##就能夠更好的查看這些基因的功能了，使用一個python腳本如下

python腳本

import pandas as pd

# 文件路徑
annotation_file = 'annotation.csv'  # CSV 格式注釋文件
gene_id_file = 'gene_id_move.txt'  # 基因 ID 文件
output_file = 'filtered_annotation.csv'  # 輸出的篩選結果文件

# 讀取 gene_id_move.txt 文件，獲取基因ID列表
with open(gene_id_file, 'r') as f:
    gene_ids = [line.strip() for line in f if line.strip()]

# 讀取 CSV 文件
df = pd.read_csv(annotation_file)

# 提取包含這些 GeneID 的行
filtered_df = df[df['GeneID'].isin(gene_ids)]

# 保存提取結果到新的 CSV 文件
filtered_df.to_csv(output_file, index=False)

print(f"提取完成，共提取 {len(filtered_df)} 條注釋信息，已保存到 {output_file}")

完成提取后，我要對功能基因進行聚類，并繪制聚類柱狀圖，并生成一個新文件用于下一步分析

import pandas as pd
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42
# 文件路徑
filtered_annotation_file = 'step10_filtered_annotation.csv'

# 讀取篩選后的注釋文件
df = pd.read_csv(filtered_annotation_file)

# 提取 GeneID 和 GO 功能注釋
df = df[['GeneID', 'GO']].dropna()  # 刪除 GO 注釋為空的行

# 將多條 GO 注釋合并為一個字符串
df['GO_text'] = df['GO'].apply(lambda x: ' '.join(x.split(';')))  # 以空格拼接不同的功能

# 向量化 GO 注釋文本
vectorizer = TfidfVectorizer(stop_words='english')
X = vectorizer.fit_transform(df['GO_text'])

# 使用 KMeans 聚類
n_clusters = 5  # 可以調整聚類數(shù)量
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
df['Cluster'] = kmeans.fit_predict(X)

# 自定義配色方案（用戶可以修改此色號）
custom_colors = ['#4b4279', '#3c6682', '#2f837f', '#45a778', '#8fc456']  # 這里可以傳入自定義顏色代碼

# 聚類結果可視化
plt.figure(figsize=(10, 6))
sns.set(style='whitegrid')  # 確保背景為白色
sns.countplot(data=df, x='Cluster', palette=custom_colors[:n_clusters])  # 使用自定義配色
plt.title('GO Annotation Clustering', fontsize=16)
plt.xlabel('Cluster Label', fontsize=14)
plt.ylabel('Number of Genes', fontsize=14)

# 導出為 PDF 格式，矢量化輸出
output_pdf = 'step10_GO_annotation_clustering.pdf'
plt.savefig(output_pdf, format='pdf', bbox_inches='tight')  # 保存為 PDF
print(f"聚類結果圖已保存為 {output_pdf}")

# 保存聚類結果 CSV 文件
df[['GeneID', 'GO_text', 'Cluster']].to_csv('step_10clustered_genes.csv', index=False)
print("聚類結果已保存到 step_10clustered_genes.csv")

image.png

import pandas as pd
from collections import Counter

# 文件路徑
csv_file = 'step_10clustered_genes.csv'
output_file = 'Cluster_key_function.txt'

# 讀取 CSV 文件
df = pd.read_csv(csv_file)

# 打開輸出文件
with open(output_file, 'w') as f:
    f.write("Cluster Key Function Report\n")
    f.write("=" * 50 + "\n\n")

    # 查看每個聚類的主要功能
    clusters = df['Cluster'].unique()
    for cluster in clusters:
        cluster_genes = df[df['Cluster'] == cluster]['GO_text']
        all_go_terms = ' '.join(cluster_genes).split()  # 將所有 GO 詞匯合并
        go_counter = Counter(all_go_terms)  # 統(tǒng)計詞頻
        most_common_go = go_counter.most_common(10)  # 統(tǒng)計前 10 個常見功能
        
        # 寫入文件
        f.write(f"\nCluster {cluster} - Top 10 GO Terms:\n")
        for term, count in most_common_go:
            f.write(f"{term}: {count}\n")
        f.write("\n" + "=" * 50 + "\n")

print(f"統(tǒng)計結果已保存到 {output_file}")

最后，繪制了一個詞云圖

import pandas as pd
from wordcloud import WordCloud
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42
# 文件路徑
csv_file = 'step_10clustered_genes.csv'
output_pdf = 'GO_terms_wordcloud.pdf'  # 輸出 PDF 文件名

try:
    # 讀取 CSV 文件
    df = pd.read_csv(csv_file)

    # 提取所有 GO 注釋文本
    all_go_text = ' '.join(df['GO_text'].dropna())

    # 創(chuàng)建詞云圖，限制最多 300 個詞
    wordcloud = WordCloud(width=800, height=400, background_color='white', max_words=300, colormap='cividis').generate(all_go_text)

    # 顯示詞云圖
    plt.figure(figsize=(8, 5))
    plt.imshow(wordcloud, interpolation='bilinear')
    plt.axis('off')
    plt.title("GO Terms Word Cloud (Top 300 Words)", fontsize=16)

    # 保存為 PDF 文件
    plt.savefig(output_pdf, format='pdf', bbox_inches='tight')
    print(f"詞云圖已保存為 {output_pdf}")

    # 顯示圖像
    plt.show()

except FileNotFoundError:
    print(f"錯誤：文件 '{csv_file}' 未找到，請檢查路徑！")
except Exception as e:
    print(f"發(fā)生錯誤：{e}")

image.png

第六部分

第五部分中，我沒有計算基因組中CDS區(qū)域和內含子區(qū)的總長度，我利用一個python腳本，從gff文件中統(tǒng)計并計算了

from collections import defaultdict

def parse_gff(gff_file):
    """
    Parses the GFF file to extract mRNA and CDS regions.

    Parameters:
    gff_file (str): Path to the GFF file.

    Returns:
    dict: Dictionary with mRNA lengths and CDS regions per gene.
    """
    genes_data = defaultdict(lambda: {'mRNA_length': 0, 'CDS_regions': []})

    with open(gff_file, 'r') as f:
        for line in f:
            if line.startswith("#"):
                continue
            columns = line.strip().split("\t")
            seq_type = columns[2]
            start, end = int(columns[3]), int(columns[4])
            attributes = columns[8]

            if seq_type == 'mRNA':
                gene_id = attributes.split("ID=")[1].split(";")[0]
                genes_data[gene_id]['mRNA_length'] = end - start + 1

            elif seq_type == 'CDS':
                parent_id = attributes.split("Parent=")[1].split(";")[0]
                genes_data[parent_id]['CDS_regions'].append((start, end))

    return genes_data


def calculate_lengths(genes_data):
    """
    Calculate total CDS and intron lengths.

    Parameters:
    genes_data (dict): Dictionary with mRNA lengths and CDS regions.

    Returns:
    tuple: Total CDS length, total intron length.
    """
    total_cds_length = 0
    total_intron_length = 0

    for gene_id, data in genes_data.items():
        cds_length = sum(end - start + 1 for start, end in data['CDS_regions'])
        total_cds_length += cds_length

        mRNA_length = data['mRNA_length']
        intron_length = mRNA_length - cds_length

        total_intron_length += max(0, intron_length)

    return total_cds_length, total_intron_length


def main():
    gff_file = "NEW_WDS.gff"
    genes_data = parse_gff(gff_file)
    cds_length, intron_length = calculate_lengths(genes_data)

    print(f"Total CDS length: {cds_length} bp")
    print(f"Total intron length: {intron_length} bp")


if __name__ == "__main__":
    main()

同樣第五部分中，老師說那個環(huán)形圖偏向性太強，看了之后好像我們的變異位點集中在內含子區(qū)，是因為沒有增加基因組各部分的比例，因此我繪制了一個雙環(huán)圖

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams

# 設置字體和輸出格式為矢量化
rcParams['pdf.fonttype'] = 42  # 使用 Type 3 字體，AI 軟件兼容性更好
rcParams['ps.fonttype'] = 42

# 讀取 Excel 文件
data_file = "step8_1_1exon.xlsx"
df = pd.read_excel(data_file)

# 檢查數(shù)據(jù)格式是否正確
#print("Data preview:")
#print(df.head())

# 計算總和
df_sum_snp = pd.DataFrame({
    "region": ["exon", "intron", "promoter_region", "other_region"],
    "SNP_count": [df["exon"].sum(), df["introns"].sum(), df["promoter region"].sum(), df["other region"].sum()]
})

# 提取區(qū)域長度數(shù)據(jù)
region_lengths = {
    "exon": df.loc[df['ID'] == 'length', 'exon'].values[0],
    "intron": df.loc[df['ID'] == 'length', 'introns'].values[0],
    "promoter_region": df.loc[df['ID'] == 'length', 'promoter region'].values[0],
    "other_region": df.loc[df['ID'] == 'length', 'other region'].values[0]
}
df_sum_length = pd.DataFrame({
    "region": ["exon", "intron", "promoter_region", "other_region"],
    "length": [region_lengths["exon"], region_lengths["intron"], region_lengths["promoter_region"], region_lengths["other_region"]]
})

# 繪制雙層環(huán)形圖函數(shù)
def draw_double_donut_chart(snp_data, length_data, colors=None, output_file="step8_1_1double_donut_chart.pdf"):
    plt.figure(figsize=(10, 10))

    # 內層環(huán)形圖 - SNP 位點占比
    plt.pie(
        snp_data["SNP_count"],
        labels=snp_data["region"],
        radius=0.7,
        autopct='%1.1f%%',
        startangle=90,
        wedgeprops={'width': 0.3, 'edgecolor': 'w'},
        colors=colors
    )

    # 外層環(huán)形圖 - 區(qū)域長度占比
    plt.pie(
        length_data["length"],
        labels=length_data["region"],
        radius=1.0,
        autopct='%1.1f%%',
        startangle=90,
        wedgeprops={'width': 0.1, 'edgecolor': 'w'},
        colors=colors
    )

    plt.title("Distribution of SNP Counts and Genomic Region Lengths", pad=20)

    # 保存為矢量化 PDF
    plt.savefig(output_file, format='pdf', bbox_inches='tight')
    plt.show()
    print(f"Chart saved as {output_file}")

# 示例顏色輸入
custom_colors = ["#fb876e", "#b4cd58", "#7dd8cf", "#7e9d31"]  # 用戶可以根

image.png

第七部分

棄掉糞便獨有的SNP位點，然后計算各種遺傳參數(shù)與血液數(shù)據(jù)再次進行比對

### 1.第一步提取出糞便獨有位點VCF文件中的染色體和坐標
 bcftools query -f'%CHROM\t%POS\n' feces_only_clean_fixed.vcf > feces_position.txt
### 2.第二步將糞便獨有位點過濾掉，生成一個新的VCF文件
bcftools view -T ^feces_position.txt -Oz -o AT_faece_MAF.filtered.vcf.gz ../../var/AT_faece_MAF.recode.vcf.gz
### 查看結果（位點數(shù)）
bcftools view -H -v snps AT_faece_MAF.filtered.vcf.gz | wc -l