过滤大型 vcf 文件答案

【问题标题】：Filtering large vcf file过滤大型 vcf 文件
【发布时间】：2022-11-05 23:13:36
【问题描述】：

我有一个具有以下格式的 VCF 文件：

#CHROM	POS	ID	REF	ALT	QUAL	FILTER.	INFO
chr1	10061	.	A	T	77.1	AC0	AC=2;AN=53780
chr1	10162	.	A	GC	81.0.	AC0;AS_VQSR	AC=1;AN=3615

我想应用几个过滤器：

仅保留 REF 和 ALT 列长度正好为 1 的行。
在第一次过滤后，我想保留那些 AC（查看 INFO）列应高于某个阈值的单元格。
最后根据 chr1 和 Pos 删除重复项，从而保留最高质量的行（QUAL 列）。
因此，如果 AC 的阈值为 2 或更多，则输出将如下所示：

#CHROM POS ID REF ALT QUAL FILTER. INFO

chr1 10061 . A T 77.1 AC0 AC=2;AN=53780

它是一个超过 845923625 行的大压缩文件。我正在考虑通过 pandas 阅读它，因为它是制表符分隔的。因此，有人可以帮助我以最有效的方式过滤此文件。谢谢！！！

【问题讨论】：

标签： python pandas

【解决方案1】：

使用以下模仿您的玩具数据框：

import pandas as pd

df = pd.DataFrame(
    {
        "#CHROM": ["chr1", "chr1", "chr2", "chr1"],
        "POS": [10061, 10162, 10163, 10061],
        "ID": [".", ".", ".", "."],
        "REF": ["A", "A", "AA", "A"],
        "ALT": ["T", "GC", "Y", "Z"],
        "QUAL": ["77.1", "81.0.", "80.0", "63.0"],
        "FILTER.": ["AC0", "AC0;AS_VQSR", "AC1", "AC2"],
        "INFO": ["AC=2;AN=53780", "AC=1;AN=3615", "AC=0;AN=3615", "AC=2;AN=3615"],
    }
)

print(df)
# Output
  #CHROM    POS ID REF ALT   QUAL      FILTER.           INFO
0   chr1  10061  .   A   T   77.1          AC0  AC=2;AN=53780
1   chr1  10162  .   A  GC  81.0.  AC0;AS_VQSR   AC=1;AN=3615
2   chr2  10163  .  AA   Y   80.0          AC1   AC=0;AN=3615
3   chr1  10061  .   A   Z   63.0          AC2   AC=2;AN=3615

这是一种方法：

df = (
    df.loc[
        (df["REF"].str.len() == 1)
        & (df["ALT"].str.len() == 1)
        & (int(df["INFO"].values[0][3]) >= 2),
        :,
    ]
    .sort_values(by="QUAL", ascending=False)
    .drop_duplicates(subset=["#CHROM", "POS"], keep="first")
)

然后：

print(df)
# Output
  #CHROM    POS ID REF ALT  QUAL FILTER.           INFO
0   chr1  10061  .   A   T  77.1     AC0  AC=2;AN=53780

【讨论】：