使用 python 绘制多个样本的 SNP 密度

Posted 2023-02-18

【中文标题】使用 python 绘制多个样本的 SNP 密度

【英文标题】：Plot SNP density for multiple samples using python

【发布时间】：2021-12-06 19:45:57

【问题描述】：

已编辑

你好

我想创建一个 python 程序，它以 FCV file、window 和 increment value 作为输入，并返回一个 plot，其中每个都有 SNP 密度所有样本（列）的窗口。 示例图片如下。

我希望采取的步骤：

建立一个 X 碱基宽的窗口并计算 该窗口中的多态性 记录多态计数和窗口的起始位置 将窗口向下移动 Y 碱基，计算窗口中的多态性数量。您将计算许多在上一个窗口中计算的相同多态性。 记录多态计数和窗口的当前起始位置 继续将窗口沿染色体向下移动 Y 碱基，计算多态性，并记录计数和位置数据，直到窗口到达染色体末端 对数据框中的所有个人执行此操作 为每个人创建（计数、位置）数据的折线图或散点图。图表应为每个人显示一条线

我可以使用 R/Bioconductor 包或 Biopython 来完成，但我需要一个基本的 python 解决方案。 请提供任何帮助！ 谢谢

这是我尝试过的：VCFfile

#!/usr/bin/env python
# libraries
import argparse
import io
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

    ## Read VCF file
# Read vcf file without headers
def read_vcf(path):
    with open(path, 'r') as f:
        lines = [l for l in f if not l.startswith('##')]
    return pd.read_csv(
        io.StringIO(''.join(lines)),
        dtype='#CHROM': str, 'POS': int, 'ID': str, 'REF': str, 'ALT': str,
               'QUAL': str, 'FILTER': str, 'INFO': str,
        sep='\t'
    ).rename(columns='#CHROM': 'CHROM')

df = read_vcf('VCFFile.vcf')

# cleaning data
## format CHROM column
df['CHROM'] = df['CHROM'].str.replace('chr0','').astype(int)

## select useful columns: all columns except not useful ones
df = df[df.columns.difference(['ID', 'INFO', 'REF', 'ALT', 'QUAL', 'FILTER', 'FORMAT'])]

# Get alleles for each sample
def get_alleles(df):
    for i in df.columns.difference(['CHROM', 'POS']):
        suffix=  str(i) + '_genotype'
        df[suffix] = df[str(i)].astype(str).str[0:3]
        #df.drop(str(i), axis=1)
        #df = df[df.columns.drop(str(i))]
# apply the function
get_alleles(df)

# remove original genotype columns
filter_col = [col for col in df if col.endswith('genotype')]
filter_col.append('CHROM')
filter_col.append('POS')

df = df[filter_col]

# replace genotypes: 1/1 by 1, else by 0
list_values = ['0/0', './.', './0', '0/.', '1/0', '0/1']
df = df.replace(to_replace =list_values, value ='NaN')
df = df.replace(to_replace ='1/1', value =1)

现在我想绘制每个样本的 SNP 密度：

# plot SNP density for each sample ==========================================
# get data for each sample
# create a function to select columns
def select_sample(col):
    x = df[['POS', str(col)]]
    #remove NaN
    x = x[x[str(col)] ==1]
    return x

sample_1 = select_sample("A_genotype")
sample_2 = select_sample("B_genotype")
sample_3 = select_sample("C_genotype")
sample_4 = select_sample("D_genotype")
sample_5 = select_sample("E_genotype")
sample_6 = select_sample("F_genotype")
sample_7 = select_sample("I_genotype")
sample_8 = select_sample("P_genotype")

我无法添加 incrementValue 以获得如下图。图 1 - 使用 1,000,000 窗口大小和 100,000 增量的多态性密度图

def plot_windowed_variant_density(pos, window_size, incrementValue=None, title, ax):

    # setup windows 
    bins = np.arange(0, pos.max(), window_size)
    print(bins)
    
    #incrementValue
    #incrementValue = ???????????
    
    # use window midpoints as x coordinate
    x = (bins[1:] + bins[:-1])/2
    
    # compute variant density in each window
    count, _ = np.histogram(sample['POS'], bins=bins)
    y= count
    # plot
    sns.despine(ax=ax, offset=10)
    ax.plot(x, y)
    ax.set_xlabel('Chromosome position (Mb)')
    ax.set_ylabel('Count')
    if title:
        ax.set_title(title)
#====================================================

fig, ax = plt.subplots(figsize=(12, 3))
# Apply the function: 
for i in [sample_1, sample_2, sample_3, sample_4, sample_5, sample_6, sample_7, sample_8]:
    plot_windowed_variant_density(i.POS, 1000000,'test', ax)

【问题讨论】：

我无法回答这个问题，因为它越来越技术化，但我认为示例数据帧的循环处理将以添加下一个循环的形式使用以下代码计算bin数并处理x轴限制，然后执行当前函数。

【参考方案1】：

如果将图形的 ax 添加到函数参数中，则可以在同一个图形上创建叠加层。

# plot SNP density ==========================================
def plot_windowed_variant_density(pos, window_size, title, ax):

    # setup windows 
    bins = np.arange(0, pos.max(), window_size)

    # use window midpoints as x coordinate
    x = (bins[1:] + bins[:-1])/2
    
    # compute variant density in each window
    count, _ = np.histogram(pos, bins=bins)

    y= count

    # plot
    sns.despine(ax=ax, offset=10)
    ax.plot(x, y)
    ax.set_xlabel('Chromosome position (Mb)')
    ax.set_ylabel('Count')
    if title:
        ax.set_title(title)
#====================================================

fig, ax = plt.subplots(figsize=(12, 3))
# Apply the function: I can use a for loop
for i in [sample_1,sample_2,sample_3]:
    plot_windowed_variant_density(i.POS, 1000000,'test', ax)
    #plot_windowed_variant_density(sample_2.POS, 1000000,'test', ax)

【讨论】：

感谢您的回复。事实上，我仍然错过了重要的一步。我想计算每个窗口的SNP个数，每次滑动10 000 bp，也就是说：取一个大小为1 Mb的窗口，计算SNP个数，步长为10 000 bp，重新计算新窗口中的 SNP 数 .... 直到染色体末端。

我知道问题是关于在同一个图表上显示多条线，我正在处理代码。我根本不是基因相关的人，所以我不明白你的评论。你能编辑一下这个问题吗？

我编辑了整个脚本。我希望很清楚。谢谢