1

我正在使用 GWAS 数据，试图选择连锁不平衡独立基因座。我的方法是将所有重要的 SNP 从最重要的 -> 最不重要的顺序排列，然后在 1KG 的列表中删除任何处于连锁不平衡 (r2 > 0.2) 且 SNP 高于它们的 SNP。

我有一个文件，其中我的 SNP 排名最高 --> 最不重要：

rs2021722
rs1117490
rs2844776
rs971570

我还有一个文件列出了 LD 中的 SNP，其中的每一个（来自 SNAP）：

SNP  Proxy  Distance    RSquared    Chromosome  Coordinate_HG18
rs2021722   rs2021722   0   1.000   chr6    30282110
rs2021722   rs885912    502 1.000   chr6    30282612
rs2021722   rs971570    1618    1.000   chr6    30280492
rs2021722   rs2844776   2304    1.000   chr6    30279806
rs2021722   rs1117490   3621    1.000   chr6    30278489
rs1117490   rs1117490   0   1.000   chr6    30278489
rs1117490   rs2517610   230 1.000   chr6    30278259
rs2844776   rs971570    686 1.000   chr6    30280492
rs2844776   rs1117490   1317    1.000   chr6    30278489
rs971570    rs2021722   1618    1.000   chr6    30282110
rs971570    rs1117490   2003    1.000   chr6    30278489

我想执行一个脚本，它将读取第一个文件中的 SNP ID，在第二个文件中找到该 SNP ID，然后读取第二个文件的“代理”列。如果第二个文件中没有任何代理 SNP 位于第一个文件的较低行号（即文件中较高的位置，具有更好的排名），我希望将该 SNP ID 写入我的输出文件。

在此示例中，我的输出文件如下所示：

rs2021722
rs117490

我在 awk 和 bash 方面有一些经验，但对两者都很陌生，不知道从哪里开始完成这项任务。非常感谢任何指针。

4

2 回答 2

2 
awk 'FNR == NR {a[$1]=++n; next}
    FNR > 1 { b[$1] = (!b[$1] || (a[$2] && a[$2]<b[$1])) ? a[$2] : b[$1] } 
    END { for(i in b) if(a[i]<=b[i]) print i }
' file1 file2

笔记：

  • 第一行，对第一个文件中的 SNP id 进行排序并将结果保存到数组“a”中
  • 第 2 行，根据第 2 列获取代理 id 的最高排名（a[$2] 的最小非空值）并将结果保存到数组“b”。（FNR>1 跳过标题行）
  • 第三行，打印满足：a[i] <= b[i] 的键
于 2014-08-08T18:38:26.023 回答
0

我不完全理解您要对输出文件做什么，但至于从第一个文件中查找第二个文件中的 SNP ID，您可以使用它。 

while read line; do
    val="$line";
    echo "val: $val";
    while read line2; do
       val2=$(echo "$line2" | awk '{print $1}');
       if [ $val == $val2 ] ; then
           echo $val "=" $val2;
           proxy=$(echo "$line2" | awk '{print $2}');
           echo "proxy = $proxy";
       fi
    done < file2.txt
done <file1.txt

基本上它遍历第一个文件，然后检查第二个文件是否存在该行，如果存在，它会抓取第二个字段（代理）。如果你能更好地解释输出文件的事情，我会试一试。我想这应该能让你继续前进。

于 2014-08-08T16:04:24.360 回答