我有一个文件包含单核苷酸多态性数据称为snp.bed
,看起来像这样:
head snp.bed
Chr17 214708483 214708484 Chr17:214708484
Chr17 214708507 214708508 Chr17:214708508
Chr17 214708573 214708574 Chr17:214708574
我还有一个名为 的文件intersect.bed
,如下所示:
head intersect.bed
Chr17 214708483 214708484 Chr17:214708484 Chr17 214706266 214710783 gene50573
Chr17 214708507 214708508 Chr17:214708508 Chr17 214706266 214710783 gene50573
Chr17 214708587 214708588 Chr17:214708580 Chr17 214706266 214710783 gene50573
我想打印出修改后的版本,snp.bed
其中包含附加到每行的额外列。如果 中 的一行与snp.bed
中 一行的前 4 列匹配intersect.bed
,那么我想打印整行,并通过与(基因名称)snp.bed
中相应行的最后一列相邻而获得额外的列。intersect.bed
或者,如果来自的行snp.bed
与来自的任何行都不匹配,intersect.bed
则与由字符串“NA”而不是基因名称组成的额外列相邻。
这是我想要的输出:
head snp.matched.bed
Chr17 214708483 214708484 Chr17:214708484 gene50573
Chr17 214708507 214708508 Chr17:214708508 gene50573
Chr17 214708573 214708574 Chr17:214708574 NA
我怎样才能做到这一点?
答案1
该解决方案假设文件的行首没有空格。与您的示例有什么不同,其中有这些空格。
awk '
{
str = $1$2$3$4;
}
FNR == NR {
arr[str] = $NF;
}
FNR != NR {
gene_name = arr[str] ? arr[str] : "NA";
print $0, gene_name;
}' intersect.bed snp.bed
输出
Chr17 214708483 214708484 Chr17:214708484 gene50573
Chr17 214708507 214708508 Chr17:214708508 gene50573
Chr17 214708573 214708574 Chr17:214708574 NA
答案2
这是使用 awk 的解决方案:
$ awk -F '\t' 'BEGIN{while(getline line<"intersect.bed") {N=split(line,a,"\t"); seen[a[1]"\t"a[2]"\t"a[3]"\t"a[4]]=a[N];}} {if(seen[$0]) {name=seen[$0];} else{name="NA"}; print $0 "\t" name}' snp.bed
Chr17 214708483 214708484 Chr17:214708484 gene50573
Chr17 214708507 214708508 Chr17:214708508 gene50573
Chr17 214708573 214708574 Chr17:214708574 NA
我假设单个制表符作为两个输入文件的分隔符。
另请注意,我将“第一个第四列”解释为“前四列”。
答案3
就我个人而言,我认为对于此类任务,最好使用“真正的”编程语言。我喜欢 Python,所以这里有一个 Python 脚本可以完成您想要的操作(它故意冗长,以便您可以理解并轻松修改它):
#!/usr/bin/env python2
# intersect.py
# Read data from the first file
snp_rows = []
with open("snp.bed", 'r') as snp_file:
for row in snp_file:
snp_rows.append(row.split())
# Read data from the second file
int_rows = []
with open("intersect.bed", 'r') as int_file:
for row in int_file:
int_rows.append(row.split())
# Compare data and compute results
results = []
for row in int_rows:
if row[:4] in snp_rows:
results.append(row[:4] + [row[-1]])
else:
results.append(row[:4] + ["NA"])
# Print the results
for row in results:
print(' '.join(row))
将其保存到文件中然后执行:
python2 intersect.py
只是为了好玩,这里是一个使用标准命令(只是grep
和cut
)的 Bash 解决方案:
while read row; do
match="$(grep -F "${row}" intersect.bed)";
if [[ -n "${match}" ]]; then
echo "${row} $(echo ${match} | cut -d' ' -f8)";
else
echo "${row} NA";
fi;
done < snp.bed
请注意,一般不建议使用 Bash 进行严格的文本处理。例如,请参阅以下帖子: