打印特定数量的核苷酸

Question

我有杂合子和纯合子病例的 vcf 统计数据，我想找到与我的 maf 文件匹配的结果。问题是 maf 文件中的参考字段不同，它排除了其他状态的核苷酸，例如如果您有一个 ref CAA 并且替代变体是 CAAAAA，那么在 maf 文件中您的 ref 将是 AAA。

所以我需要一个代码来更改我的文件中的 ref 字段和 alt 统计信息（可以添加单独的列 ref2 和 alt2）

这是我的文件的一个片段：

CHR POS ID REF ALT chr11 71579744 rs71049992 A ACAGCAGCTGGACTGGGAGCAGCAGGACCTG（插入案例）

chr11 124880551 rs71859853 CCGGAGT C（删除案例）

我想我应该首先计算第 4 列和第 5 列中的核苷酸数量。然后如果第 4 列中的数字大于 5（意味着删除），那么在我的 ref2 中，该位置将从下一个不同于备选核苷酸的核苷酸开始.

对于插入，我将更改一个 alt 站点并跳过参考核苷酸

因此，我想要这个：

CHR POS ID REF ALT REF2 ALT2 chr11 71579744 rs71049992 A ACAGCAGCTGGACTGGGAGCAGCAGGACCTG A CAGCAGCTGGACTGGGAGCAGCAGGACCTG

chr11 124880551 rs71859853 CCGGAGT C CGGAGT C

非常感谢您！

Answer 1

I think I should first count numbers of nucleotides in column4 and 5…

有了awk，您可以使用length函数来计算核苷酸的数量:

awk 'NR==1 {print [=10=]" REF2 ALT2"}       # assuming first line has column headers
     NR>1  {if (length()<length()) print [=10=], , gensub(, "", 1, )
            else                       print [=10=], gensub(, "", 1, ), }' file

打印特定数量的核苷酸

Printing a specific number of nucleotides

reference

vcf-vcard

maf