Подсчитать общее уникальное количество записей на основе столбцов из нескольких файлов VCF

У меня около 200 файлов с длинными строками заголовков, которые начинаются с символа «#», а затем записываются с 4 столбцами, например:

      file_1.vcf
##some_comments that span many lines
##some_comments that span many lines
#CRHOM    POS    REF    ALT
chr1   111    A       G
chr2   222    G       T
chrY   99999  A       C

file_2.vcf
##some_comments that span many lines
##some_comments that span many lines
#CRHOM    POS    REF    ALT
chr2   444    G       T
chr2   222    G       T
chrY   7473  C       A

Для этого простого примера всего с 2 vcf ожидаемое общее количество записей составляет 5 записей.

Каждый из файлов содержит около миллиона записей, и я пытаюсь подсчитать, сколько их общего количества записей с одним условием, что повторяющаяся запись, появляющаяся в любом другом файле vcf, должна учитываться только один раз. Как я могу это сделать, я пытаюсь сделать это на Python, но это застряло очень долго:

          import os
    from pathlib import Path
    import vcf

    path_gdc_data = Path(os.getcwd()+"/all_vcfs/")

    non_dup = []
    count += 1
    for i, vcf_file in enumerate(path_gdc_data.rglob("*.vcf")):
        vcf_reader = vcf.Reader(filename=str(vcf_file))
        for rec in vcf_reader:
            if (rec.CHROM,rec.POS, rec.REF, rec.ALT) not in indels_coord_uniq:
                    non_dup_records.append((rec.CHROM, rec.POS))
                    count += 1
    print(f"Total num of non dup in all vcf are {count}")

Примечание 1: я включаю теги pandas и Dataframe в свой вопрос, чтобы увидеть, будут ли pandas более эффективным решением для этого.

Примечание 2: в каждом отдельном файле не может быть дублированных записей, и каждый файл сортируется.

Пример файла vcf:

      ##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS     ID        REF ALT    QUAL FILTER INFO                              FORMAT      NA00001        NA00002        NA00003
20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
20     17330   .         T      A       3    q10    NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3   0/0:41:3
20     1110696 rs6040355 A      G,T     67   PASS   NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2   2/2:35:4
20     1230237 .         T      .       47   PASS   NS=3;DP=13;AA=T                   GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
20     1234567 microsat1 GTCT   G,GTACT 50   PASS   NS=3;DP=9;AA=G                    GT:GQ:DP    0/1:35:4       0/2:17:2       1/1:40:3

Спасибо заранее.

Источник

2 ответа

Если файлы vcf отсортированы и в одном файле нет уникальных пользователей, вы можете использовать heapq.merge чтобы объединить все файлы вместе, а затем использовать itertools.groupbyчтобы сгруппировать их (я полагаю, мы хотим подсчитать уникальных посетителей во всех столбцах). Затем подсчитайте количество уникальных групп. Пример

      import vcf
from heapq import merge
from itertools import groupby

filelist = ["file_1.vcf", "file_2.vcf"]  # <--- create the list from `glob()` for example
vcfs = [vcf.Reader(filename=f) for f in filelist]

uniq = sum(1 for _ in groupby(merge(*vcfs)))

print("Total number of unique records:", uniq)

Печать (в моем примере):

      Total number of unique records: 7

РЕДАКТИРОВАТЬ: для подсчета уникальных посетителей только в столбцах CHROM, POS, REF, ALT:

      import vcf
from heapq import merge
from itertools import groupby

filelist = ["file_1.vcf", "file_2.vcf"]
vcfs = [
    ((row.CHROM, row.POS, row.REF, row.ALT) for row in vcf.Reader(filename=f))
    for f in filelist
]

uniq = sum(1 for _ in groupby(merge(*vcfs)))

print("Total number of unique records:", uniq)
Источник

Другой метод с использованием Pandas:

      import pandas as pd
import io
from pathlib import Path

COLS = ["CRHOM", "POS", "REF", "ALT"]

def read_vcf(filepath_or_buffer, usecols=None):
    with open(filepath_or_buffer) as vcf:
        while True:
            line = vcf.readline()
            if not line.startswith("##"):
                names = line[1:].split()
                break
        return pd.read_table(vcf, sep="\t", delim_whitespace=True,
                             usecols=usecols, header=None, names=headers)


# Get vcf file list and load the first to a dataframe
path_gdc_data = Path().cwd() / "all_vcfs"
vcf_files = path_gdc_data.glob("*.vcf")
df = read_vcf(next(vcf_files), usecols=COLS)

# Compare and merge other dataframes
for vcf in vcf_files:
    df1 = read_vcf(vcf, usecols=COLS)
    df = df.merge(df1, on=COLS, how="outer")

print(f"Total num of non dup in all vcf are {len(df)}")
Источник
Другие вопросы по тегам