如何将基因组 SNP 转换为二进制格式？答案

【问题标题】：How to convert genomic SNPs to a binary format?如何将基因组 SNP 转换为二进制格式？
【发布时间】：2021-01-20 03:14:48
【问题描述】：

我正在努力弄清楚如何转换如下所示的 (SNP) 文件：

pos,sample1,sample2,sample3,sample4,sample5,sample6,sample7,sample8,sample9,sample10,sample11,sample12,sample13,sample14
79107,C,C,C,C,C,C,C,C,C,C,C,C,G,G
79115,C,C,C,C,C,C,C,A,C,C,T,C,C,C
79116,A,A,A,A,A,A,A,A,A,A,A,A,A,A
79124,C,C,T,T,C,C,C,C,C,C,C,C,C,C
79128,G,G,G,G,A,G,G,G,G,C,G,G,G,G

转成这样的二进制格式：

pos,sample1,sample2,sample3,sample4,sample5,sample6,sample7,sample8,sample9,sample10,sample11,sample12,sample13,sample14
79107,0,0,0,0,0,0,0,0,0,0,0,0,1,1
79115,0,0,0,0,0,0,0,1,0,0,1,0,0,0
79116,0,0,0,0,0,0,0,0,0,0,0,0,0,0
79124,0,0,1,1,0,0,0,0,0,0,0,0,0,0
79128,0,0,0,0,1,0,0,0,0,1,0,0,0,0

有谁知道 R 中的一个包，或者命令行工具，awk 代码，可以用来做这个吗？

以下解决方案的一个问题是都这样做：

1097023,A,A,G,G,G,G,G,G,G,G,G,G,A,A
1097027,C,C,C,C,C,C,C,C,C,C,C,C,C,C
4363243,C,A,A,A,A,A,A,A,A,A,A,A,A,A 
4363270,T,T,T,T,T,T,T,T,T,T,T,T,T,T 
4363275,A,G,G,G,G,G,G,G,G,G,G,G,G,G

1097023,0,0,1,1,1,1,1,1,1,1,1,1,0,0
1097027,0,0,0,0,0,0,0,0,0,0,0,0,0,0
4363243,0,1,1,1,1,1,1,1,1,1,1,1,1,1 
4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
4363275,0,1,1,1,1,1,1,1,1,1,1,1,1,1

虽然应该是：

1097023,1,1,0,0,0,0,0,0,0,0,0,0,1,1
1097027,0,0,0,0,0,0,0,0,0,0,0,0,0,0
4363243,1,0,0,0,0,0,0,0,0,0,0,0,0,0 
4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 
4363275,1,0,0,0,0,0,0,0,0,0,0,0,0,0

所以我认为解决方案应该包括如果大多数样本'即> 7个样本是一样的，这个变成0，不符合多数的变成1。

【问题讨论】：

这里的逻辑是什么？ C 在倒数第二行变成0，在最后一行变成1。我们需要有关您想要的目标格式的更多信息。
DNA 由 A、C、T 或 G 组成，因此您可以有 4 个选项。通常，您会在每个样本中看到相同的结果，因此所有 A、所有 C、所有 T 或所有 G。但是，如果您有“突变 = 从例如 C 到 G 的转换”，这称为 SNP。对于下游分析，我需要一个二进制格式，其中所有内容都为 0，除了一个或多个样本与给定位置的其余样本不匹配的地方（第一列“pos”）。
我们所要做的就是您在问题中提供的信息。发布答案时，您的问题中没有任何内容可以表明您现在添加了I think the solution should include that if the majority of samples 'i.e. > 7 samples' is the same, this becomes a 0, the ones that don't match the majority become 1.。您可能希望将这个问题恢复原状，接受对您提出的问题的回答，然后根据您的新要求和更具代表性的示例提出后续问题。

标签： python r awk bioinformatics

【解决方案1】：

这是一个简单的awk（标准Linux awk/gawk）脚本：

script.awk

BEGIN {FS = ","} # set field seperator to ","
NR>1{     # every line except the first line
  zeroIndc=$2; # identify the zero indicator and save the variable
  for (i = 2; i <= NF; i++) { # for each DNA letter
    if ($i == zeroIndc) { # if DNA letter match zero indicator
      $i = 0; # set DNA letter to 0
    } else { # if DNA letter not match zero indicator
      $i = 1; # set DNA letter to 1
    }
  }
  print; # print new line at end 
}

运行脚本`script.awk`

awk -f scirpt.awk input.txt

运行单行 awk 脚本

awk 'BEGIN{FS=","}NR>1{z=$2;for(i=2;i<=NF;i++)$i=($i==z)?0:1;print;}' input.txt

【讨论】：

这很完美，我在一个大型数据集上尝试过，它可以工作。非常感谢！
您也必须设置OFS=","，否则它会将每个逗号替换为空白。

【解决方案2】：

使用基数 R，您可以将每个 sample 列与第一个样本列进行比较，并在值不匹配的地方返回 1。

cols <- grep('sample', names(df))
df[cols] <- +(df$sample1 != df[cols])
df

#    pos sample1 sample2 sample3 sample4 sample5 sample6 sample7 sample8 sample9
#1 79107       0       0       0       0       0       0       0       0       0
#2 79115       0       0       0       0       0       0       0       1       0
#3 79116       0       0       0       0       0       0       0       0       0
#4 79124       0       0       1       1       0       0       0       0       0
#5 79128       0       0       0       0       1       0       0       0       0

#  sample10 sample11 sample12 sample13 sample14
#1        0        0        0        1        1
#2        0        1        0        0        0
#3        0        0        0        0        0
#4        0        0        0        0        0
#5        1        0        0        0        0

虽然上面的方法在处理大型数据集时会更高效，但这里是 dplyr 库的替代方法。

library(dplyr)
df <- df %>% mutate(across(contains('sample'), ~+(sample1 != .)))

数据

df <- structure(list(pos = c(79107L, 79115L, 79116L, 79124L, 79128L
), sample1 = c("C", "C", "A", "C", "G"), sample2 = c("C", "C", 
"A", "C", "G"), sample3 = c("C", "C", "A", "T", "G"), sample4 = c("C", 
"C", "A", "T", "G"), sample5 = c("C", "C", "A", "C", "A"), sample6 = c("C", 
"C", "A", "C", "G"), sample7 = c("C", "C", "A", "C", "G"), sample8 = c("C", 
"A", "A", "C", "G"), sample9 = c("C", "C", "A", "C", "G"), sample10 = c("C", 
"C", "A", "C", "C"), sample11 = c("C", "T", "A", "C", "G"), sample12 = c("C", 
"C", "A", "C", "G"), sample13 = c("G", "C", "A", "C", "G"), sample14 = c("G", 
"C", "A", "C", "G")), class = "data.frame", row.names = c(NA, -5L))

【讨论】：

我试过你的两个solotions。它们很有用，但是我遇到与发布的 awk 解决方案相同的问题：4363243,C,A,A,A,A,A,A,A,A,A,A,A,A,A 4363270,T,T,T,T,T,T,T,T,T,T,T,T,T,T 4363275,A,G,G,G,G,G,G,G,G,G,G,G,G,G 变为 4363243,0,1,1,1,1,1,1,1,1,1,1,1,1,1 4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 4363275,0,1,1,1,1,1,1,1,1,1,1,1,1,1 而不是 4363243,1,0,0,0,0,0,0,0,0,0,0,0,0,0 4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 4363275,1,0,0,0,0,0,0,0,0,0,0,0,0,0 我当然可以手动删除它们，但我想知道是否有办法解决这个问题
@XelaVi 抱歉，由于代码格式不正确，因此不清楚您想从评论中显示什么。您可以更新您的帖子以包含任何其他详细信息。

【解决方案3】：

$ awk -F, 'NR>1{gsub($2,0); gsub(/[ACTG]/,1)} 1' file
pos,sample1,sample2,sample3,sample4,sample5,sample6,sample7,sample8,sample9,sample10,sample11,sample12,sample13,sample14
79107,0,0,0,0,0,0,0,0,0,0,0,0,1,1
79115,0,0,0,0,0,0,0,1,0,0,1,0,0,0
79116,0,0,0,0,0,0,0,0,0,0,0,0,0,0
79124,0,0,1,1,0,0,0,0,0,0,0,0,0,0
79128,0,0,0,0,1,0,0,0,0,1,0,0,0,0

【讨论】：

这个和其他 awk 解决方案的唯一问题是：4363243,C,A,A,A,A,A,A,A,A,A,A,A,A,A 4363270,T,T,T,T,T,T,T,T,T,T,T,T,T,T 4363275,A,G,G,G,G,G,G,G,G,G,G,G,G,G 变为 4363243,0,1,1,1,1,1,1,1,1,1,1,1,1,1 4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 4363275,0,1,1,1,1,1,1,1,1,1,1,1,1,1 而不是 4363243,1,0,0,0,0,0,0,0,0,0,0,0,0,0 4363270,0,0,0,0,0,0,0,0,0,0,0,0,0,0 4363275,1,0,0,0,0,0,0,0,0,0,0,0,0,0 但我可以手动删除这些情况;)

【解决方案4】：

识别参考等位基因，在这种情况下，参考基于最常见的等位基因，使用data from Ronak's answer，我们有 5 个 SNP REF 等位基因：

REF <- apply(df[, -1], 1, function(i) names(which.max(table(i)))) 
# [1] "C" "C" "A" "C" "G"

现在，将每个样本的列与 REF 进行比较，然后通过将逻辑 TRUE/FALSE 乘以 1 转换为整数。并将其分配回原始 dataframe：

df[, -1] <- (df[, -1] != REF) * 1
df
#     pos sample1 sample2 sample3 sample4 sample5 sample6 sample7 sample8 sample9 sample10 sample11 sample12 sample13 sample14
# 1 79107       0       0       0       0       0       0       0       0       0        0        0        0        1        1
# 2 79115       0       0       0       0       0       0       0       1       0        0        1        0        0        0
# 3 79116       0       0       0       0       0       0       0       0       0        0        0        0        0        0
# 4 79124       0       0       1       1       0       0       0       0       0        0        0        0        0        0
# 5 79128       0       0       0       0       1       0       0       0       0        1        0        0        0        0

【讨论】：

【解决方案5】：

假设您想将“reference”等位基因编码为0，其中“reference”是最常见的，这是使用 Python 3 的解决方案，它的 pandas 数据处理库（使用 @987654324 安装@) 和Counter specialized dictionary。

从包含您的示例数据的文件 snps_letters.csv 开始：

$ cat snps_letters.csv
pos,sample1,sample2,sample3,sample4,sample5,sample6,sample7,sample8,sample9,sample10,sample11,sample12,sample13,sample14
79107,C,C,C,C,C,C,C,C,C,C,C,C,G,G
79115,C,C,C,C,C,C,C,A,C,C,T,C,C,C
79116,A,A,A,A,A,A,A,A,A,A,A,A,A,A
79124,C,C,T,T,C,C,C,C,C,C,C,C,C,C
79128,G,G,G,G,A,G,G,G,G,C,G,G,G,G

我们将它打开，转换并保存在脚本binarize.py：

#!/usr/bin/env python3

from collections import Counter
import pandas as pd

# Load the comma-separated data in a pandas DataFrame object:
snps = pd.read_csv("snps_letters.csv", index_col="pos")

def binarize(column):
    # Extracting the most common element as ref
    # (using list and tuple unpacking,
    # because the most_common method returns
    # a list of (element, count) pairs):
    [(ref, _)] = Counter(column).most_common(1)
    # We can now define an auxiliary function
    # to recode individual elements in the column:
    def to_bin(elem):
        # int(True) -> 1, int(False) -> 0
        return 1 - int(elem == ref)
    # Transform the column by applying the recoding function to its elements
    # (see https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.apply.html#pandas.Series.apply):
    return column.apply(to_bin)

# Now apply our column-transforming function to the columns (axis=0)
# (see https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.apply.html):
binarized = snps.apply(binarize, axis=0)
# Write the result to a file:
binarized.to_csv("snps_binary.csv")

在shell中执行脚本：

$ ./binarize.py

检查保存的结果：

$ cat snps_binary.csv 
pos,sample1,sample2,sample3,sample4,sample5,sample6,sample7,sample8,sample9,sample10,sample11,sample12,sample13,sample14
79107,0,0,0,0,0,0,0,0,0,0,0,0,0,0
79115,0,0,0,0,0,0,0,1,0,0,1,0,1,1
79116,1,1,1,1,1,1,1,1,1,1,1,1,1,1
79124,0,0,1,1,0,0,0,0,0,0,0,0,1,1
79128,1,1,1,1,1,1,1,1,1,0,1,1,0,0

为了有一个更可重用的脚本，你可以将输入和输出文件路径作为参数传递给命令行，它们可以在 python 代码中作为sys.argv[1] 和sys.argv[2] 访问（前提是@987654334 @先）。

在命令行上获取输入和输出文件的版本，并将参考核苷酸编码为 1 而不是 0：

#!/usr/bin/env python3

import sys
from collections import Counter
import pandas as pd

# Load the comma-separated data in a pandas DataFrame object:
snps = pd.read_csv(sys.argv[1], index_col="pos")

def binarize(column):
    # Extracting the most common element as ref
    # (using list and tuple unpacking,
    # because the most_common method returns
    # a list of (element, count) pairs):
    [(ref, _)] = Counter(column).most_common(1)
    # We can now define an auxiliary function
    # to recode individual elements in the column:
    def to_bin(elem):
        # int(True) -> 1, int(False) -> 0
        return int(elem == ref)
    # Transform the column by applying the recoding function to its elements
    # (see https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.apply.html#pandas.Series.apply):
    return column.apply(to_bin)

# Now apply our column-transforming function to the columns (axis=0)
# (see https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.apply.html):
binarized = snps.apply(binarize, axis=0)
# Write the result to a file:
binarized.to_csv(sys.argv[2])

它的用法如下：

./binarize.py snps_letters.csv snps_binary.csv

【讨论】：

嘿，我尝试运行你的脚本，得到的结果和你一样，但这并不完全正确。 79128,1,1,1,1,1,1,1,1,1,0,1,1,0,0 应该有 0,0,0,0,0,0,0,0,0,1 ,0,1,1
这发生在我的整个数据集中。我尝试设置 most_common(0)，但没有帮助。
@XelaVi 如果您希望将最常见的编码为 1 而不是 0，请将 return 1 - int(elem == ref) 更改为 return int(elem == ref)。 Counter 的 most_common 方法中的 1 指定我们只需要最常见的元素。如果我们输入2，这意味着我们想要最常见的 2 个。我不知道如果我们输入 0 会发生什么...
@XelaVi 我在末尾添加了一个新版本，用于执行反向编码。但我认为与您的期望可能仍然存在差异。对于某些样本来说，最常见的情况下存在 ex-aequo 使得选择有些武断。

script.awk

运行脚本script.awk

运行单行 awk 脚本

运行脚本`script.awk`