VCF 파일에서 missing allele 분석하는 방법

VCF 파일에서 missing allele 분석하는 방법

 

VCF 파일에서는 missing allele이 . 또는 ./. 등의 형식으로 표시됩니다.

이 정보를 추출하여 분석하는 방법을 단계별로 설명하겠습니다.

 

1. VCF 파일 읽기 및 데이터프레임 생성

import pandas as pd
import glob

# VCF 파일 경로
vcf_file = 'path/to/your/file.vcf'

# VCF 파일 읽기 함수
def read_vcf(file):
    with open(file, 'r') as f:
        lines = f.readlines()

    header_line = [line for line in lines if line.startswith('#') and not line.startswith('##')][0]
    header = header_line.strip().split('\t')
    data_lines = [line for line in lines if not line.startswith('#')]

    data = [line.strip().split('\t') for line in data_lines]
    df = pd.DataFrame(data, columns=header)
    return df

# VCF 파일 읽기
df = read_vcf(vcf_file)
print(df.head())

 

2. GT (Genotype) 필드 추출 및 결측치 분석

# FORMAT 필드 파싱 함수
def parse_format(format_str, sample_str):
    format_keys = format_str.split(':')
    sample_values = sample_str.split(':')
    return dict(zip(format_keys, sample_values))

# 모든 샘플의 GT 필드를 추출
sample_columns = df.columns[9:]  # 샘플 열 (9번째 열 이후)
gt_data = []

for index, row in df.iterrows():
    format_str = row['FORMAT']
    gt_dicts = [parse_format(format_str, row[col]) for col in sample_columns]
    gt_data.append(gt_dicts)

# GT 필드를 데이터프레임으로 변환
gt_df = pd.DataFrame(gt_data, columns=sample_columns)
gt_df = gt_df.applymap(lambda x: x['GT'] if 'GT' in x else None)

# GT 데이터프레임 확인
print(gt_df.head())

 

3. Missing allele 분석

GT 데이터프레임에서 missing allele(. 또는 ./.)의 개수를 계산합니다.

# Missing allele을 `.` 또는 `./.`로 간주하고, 각 샘플별 결측치 개수 계산
missing_alleles = gt_df.applymap(lambda x: x in ['.', './.'])
missing_counts = missing_alleles.sum()

# 각 샘플별 결측치 개수 출력
print(missing_counts)

 

4. Missing allele 비율 계산

각 샘플에서 missing allele의 비율을 계산합니다.

# 각 샘플의 전체 유전자형 개수
total_counts = gt_df.count()

# 결측치 비율 계산
missing_ratio = missing_counts / total_counts

# 결측치 비율 출력
print(missing_ratio)

 

5. 시각화

Matplotlib를 사용하여 각 샘플에서의 missing allele 비율을 시각화합니다.

import matplotlib.pyplot as plt

# 결측치 비율 시각화
plt.figure(figsize=(12, 6))
missing_ratio.plot(kind='bar')
plt.title('Missing Allele Ratio per Sample')
plt.xlabel('Sample')
plt.ylabel('Missing Allele Ratio')
plt.show()

 

통합 코드

import pandas as pd
import matplotlib.pyplot as plt
import glob

# VCF 파일 경로
vcf_file = 'path/to/your/file.vcf'

# VCF 파일 읽기 함수
def read_vcf(file):
    with open(file, 'r') as f:
        lines = f.readlines()

    header_line = [line for line in lines if line.startswith('#') and not line.startswith('##')][0]
    header = header_line.strip().split('\t')
    data_lines = [line for line in lines if not line.startswith('#')]

    data = [line.strip().split('\t') for line in data_lines]
    df = pd.DataFrame(data, columns=header)
    return df

# FORMAT 필드 파싱 함수
def parse_format(format_str, sample_str):
    format_keys = format_str.split(':')
    sample_values = sample_str.split(':')
    return dict(zip(format_keys, sample_values))

# VCF 파일 읽기
df = read_vcf(vcf_file)

# 모든 샘플의 GT 필드를 추출
sample_columns = df.columns[9:]  # 샘플 열 (9번째 열 이후)
gt_data = []

for index, row in df.iterrows():
    format_str = row['FORMAT']
    gt_dicts = [parse_format(format_str, row[col]) for col in sample_columns]
    gt_data.append(gt_dicts)

# GT 필드를 데이터프레임으로 변환
gt_df = pd.DataFrame(gt_data, columns=sample_columns)
gt_df = gt_df.applymap(lambda x: x['GT'] if 'GT' in x else None)

# Missing allele을 `.` 또는 `./.`로 간주하고, 각 샘플별 결측치 개수 계산
missing_alleles = gt_df.applymap(lambda x: x in ['.', './.'])
missing_counts = missing_alleles.sum()

# 각 샘플의 전체 유전자형 개수
total_counts = gt_df.count()

# 결측치 비율 계산
missing_ratio = missing_counts / total_counts

# 결측치 비율 출력
print(missing_ratio)

# 결측치 비율 시각화
plt.figure(figsize=(12, 6))
missing_ratio.plot(kind='bar')
plt.title('Missing Allele Ratio per Sample')
plt.xlabel('Sample')
plt.ylabel('Missing Allele Ratio')
plt.show()

 

https://alimanfoo.github.io/2017/06/14/read-vcf.html

Extracting data from VCF files