linux anaconda + bcftools使用

本来就是想在ubuntu里用一下虚拟环境 

结果miniconda一直出问题

主要就是装好了创建不了带python的虚拟环境

报错：

Verifying transaction: failed

Proceed ([y]/n)? y


Downloading and Extracting Packages:

Preparing transaction: done
Verifying transaction: failed

# >>>>>>>>>>>>>>>>>>>>>> ERROR REPORT <<<<<<<<<<<<<<<<<<<<<<

                                    ......
    OSError: [Errno 40] Too many levels of symbolic links: '/mnt/d/linuxConda/pkgs/ncurses-6.4-h6a678d5_0/share/terminfo/n/ncr260vt300wpp'

`$ /mnt/d/linuxConda/bin/conda create -n zky python==3.9.12 -c defaults`

而且 运行安装 速度奇慢

各种方法都试了 现在也没找到原因 

最后的解决方法是 

不指定虚拟机里conda的安装路径到其他盘了 而是直接安装在win11 WSL默认目录下 并安装conda而不是miniconda

————————————————————————————————

bcftools使用

1.表头问题

(zky) zky@zky:/mnt/e/pyProjects/zc/snpEvaluate$ bcftools stats test.vcf > stats.txt
[W::bcf_hdr_check_sanity] PL should be declared as Number=G
[W::vcf_parse] Contig '1A' is not defined in the header. (Quick workaround: index the file with tabix.)

才知道原来vcf的表头是有用的
修改了一下 添加了这几行

##contig=<ID=1A>
##contig=<ID=1B>
##contig=<ID=1D>
##contig=<ID=2A>
##contig=<ID=2B>
##contig=<ID=2D>
##contig=<ID=3A>
##contig=<ID=3B>
##contig=<ID=3D>
##contig=<ID=4A>
##contig=<ID=4B>
##contig=<ID=4D>
##contig=<ID=5A>
##contig=<ID=5B>
##contig=<ID=5D>
##contig=<ID=6A>
##contig=<ID=6B>
##contig=<ID=6D>
##contig=<ID=7A>
##contig=<ID=7B>
##contig=<ID=7D>

2.生成含有包括GT（基因型）、AD（等位基因深度）和DP（读取深度）信息的vcf文件

安装bcftool之后
 

sudo apt install bcftools

linux~ubuntu：

(bcftools view -h test.vcf; bcftools +fill-tags test.vcf -Ov | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t%INFO\tFORMAT[\t%GT]\n') > output.vcf


解释：

——————————————————————————————————————————————————————————
(...): 括号用于将多个命令组合在一起，作为一个整体执行。
bcftools view -h test.vcf:

bcftools view: 用于查看和过滤VCF/BCF文件。
-h: 只输出头部信息（header）。
test.vcf: 输入的VCF文件名。


;: 用于分隔两个命令。
bcftools +fill-tags test.vcf -Ov:

+fill-tags: bcftools的一个插件，用于计算和添加缺失的标签（如AF）。
test.vcf: 输入的VCF文件名。
-Ov: 输出格式为未压缩的VCF（O表示输出，v表示VCF格式）。


|: 管道符，将前一个命令的输出作为后一个命令的输入。
bcftools query -f '...':

query: 用于提取VCF文件中的特定字段。
-f: 指定输出格式。


格式字符串 '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t%INFO\tFORMAT[\t%GT]\n':

%CHROM: 染色体
%POS: 位置
%ID: 变异ID
%REF: 参考等位基因
%ALT: 替代等位基因
%QUAL: 质量分数
%FILTER: 过滤器状态
%INFO: INFO字段内容
FORMAT: FORMAT字段
[\t%GT]: 对每个样本，输出制表符和基因型
\n: 换行符


> output.vcf: 将结果重定向到名为output.vcf的文件。

总的来说，这个命令做了以下几件事：

提取原始VCF文件的头部信息。
使用fill-tags插件添加缺失的标签。
查询并格式化输出指定的字段，包括每个样本的基因型。
将头部信息和格式化的数据行组合在一起，输出到一个新的VCF文件。

3.计算每个位点的MAF、遗传多样性和PIC值，然后计算每条染色体和基因组的平均值。

PIC的正确公式是：
对于双等位基因座：
PIC = 1 - Σ(pi^2) - Σ(Σ(2pi^2pj^2))
其中：

pi 是第i个等位基因的频率
pj 是第j个等位基因的频率 (j ≠ i)

对于仅有两个等位基因的情况（这通常是SNP的情况），公式可以简化为：
PIC = 1 - p^2 - q^2 - 2p^2q^2
其中：

p 是主要等位基因的频率
q 是次要等位基因的频率 (q = 1 - p)
————————————————————————————————————————————————————————————————————————————
预期杂合度（Expected Heterozygosity）【基因多样性（Gene Diversity）】：
公式是：
H = 2p(1-p)
其中：

H 是遗传多样性
p 是一个等位基因的频率
1-p 就是另一个等位基因的频率（对于双等位基因座）

在这个函数中，我们使用 MAF（Minor Allele Frequency，次要等位基因频率）作为输入。因此：

maf 代表次要等位基因的频率
1 - maf 就是主要等位基因的频率

解释：

这个公式计算的是在随机交配群体中，从该位点随机抽取两个等位基因时，这两个等位基因不同的概率。
2 * maf * (1 - maf) 实际上计算了两种可能的异质结合子的概率之和：

maf * (1-maf)：次要等位基因和主要等位基因的组合
(1-maf) * maf：主要等位基因和次要等位基因的组合


这两种情况的概率是相同的，所以我们将其乘以2。
遗传多样性的值范围从0到0.5：

当 MAF = 0 或 1 时（即位点是单态的），H = 0
当 MAF = 0.5 时（即两个等位基因频率相等），H 达到最大值0.5


较高的遗传多样性值表示群体中该位点的变异程度较高。

用于计算的python代码：

import pandas as pd


def calculate_genetic_diversity(maf):
    return 2 * maf * (1 - maf)


def calculate_pic(maf):
    p = 1 - maf
    q = maf
    return 1 - (p**2) - (q**2) - (2 * p**2 * q**2)


def process_vcf(vcf_file):
    # 读取VCF文件，只选择需要的列
    df = pd.read_csv(vcf_file, sep='\t', comment='#',
                     usecols=['CHROM', 'POS', 'INFO'],
                     dtype={'CHROM': str, 'POS': int, 'INFO': str})

    # 重命名'#CHROM'列为'CHROM'
    df = df.rename(columns={'CHROM': 'CHROM'})

    # 从INFO列提取所需信息
    df['MAF'] = df['INFO'].str.extract(r'MAF=([\d\.]+)').astype(float)
    df['NS'] = df['INFO'].str.extract(r'NS=(\d+)').astype(int)
    df['AN'] = df['INFO'].str.extract(r'AN=(\d+)').astype(int)

    # 计算Genetic Diversity和PIC
    df['Genetic_Diversity'] = calculate_genetic_diversity(df['MAF'])
    df['PIC'] = calculate_pic(df['MAF'])

    return df[['CHROM', 'POS', 'MAF', 'NS', 'AN', 'Genetic_Diversity', 'PIC']]


def summarize_by_genome(df):
    # 定义染色体组
    a_genome = ['1A', '2A', '3A', '4A', '5A', '6A', '7A']
    b_genome = ['1B', '2B', '3B', '4B', '5B', '6B', '7B']
    d_genome = ['1D', '2D', '3D', '4D', '5D', '6D', '7D']
    all_chroms = a_genome + b_genome + d_genome

    # 计算每条染色体的平均值
    chrom_means = df.groupby('CHROM')[['MAF', 'Genetic_Diversity', 'PIC']].mean()

    # 计算基因组平均值
    genome_means = pd.DataFrame({
        'A genome': df[df['CHROM'].isin(a_genome)][['MAF', 'Genetic_Diversity', 'PIC']].mean(),
        'B genome': df[df['CHROM'].isin(b_genome)][['MAF', 'Genetic_Diversity', 'PIC']].mean(),
        'D genome': df[df['CHROM'].isin(d_genome)][['MAF', 'Genetic_Diversity', 'PIC']].mean(),
    }).T

    # 计算总体平均值
    overall_mean = df[['MAF', 'Genetic_Diversity', 'PIC']].mean().to_frame().T
    overall_mean.index = ['Mean']

    # 合并所有结果
    results = pd.concat([chrom_means.loc[all_chroms], genome_means, overall_mean])

    return results


# 主执行
vcf_file = 'bcf_output1.vcf'
df = process_vcf(vcf_file)
summary = summarize_by_genome(df).round(3)
print(summary)

# 保存结果到CSV
summary.to_csv('genome_summary.csv', sep="\t", index_label="chrom")