特殊.bed到.fa转换(基因组坐标/.fa序列)

Question

我的目标是从基因组坐标(origin.bed)创建一个定制的蛋白质序列参考文件(origin.bed)。

(origin.bed；与染色体，起始，结束，TranscriptID，链，GeneID)

chr1    109202569       109202584       ENST00000370031.1_uORF_0        -       ENSG00000162639.11
chr1    109203584       109203617       ENST00000370031.1_uORF_0        -       ENSG00000162639.11
chr11   102188276       102188302       ENST00000263464.3_uORF_0        +       ENSG00000023445.9
chr11   10830291        10830306        ENST00000530211.1_uORF_1        -       ENSG00000110321.11
chr11   10830400        10830451        ENST00000530211.1_uORF_0        -       ENSG00000110321.11

期望的两步输出:基因组核苷酸序列文件(nucleotide.fa；0)和蛋白质序列文件(protein.fa)。

(插图；序列不适合于片段中的坐标)

>chr1:109202569-109203617(ENST00000370031.1_uORF_0)
TCACCCGGAACAGATTTTTTCCTCCTTGGATGACTCCGTCATGGACAC
>chr11:102188276-102188302(ENST00000263464.3_uORF_0)
nucleotide sequence xyz
>chr11:10830291-10830306(ENST00000263464.3_uORF_0)
nucleotide sequence xyz
>chr11:10830400-10830451(ENST00000263464.3_uORF_0)
nucleotide sequence xyz

>chr1:109202569-109203617(ENST00000370031.1_uORF_0)
SPGTDFFLLGULRHGH
>chr11:102188276-102188302(ENST00000263464.3_uORF_0)
protein sequence xyz
...

重要的一点是，具有相同TranscriptID的条目应该在蛋白质序列翻译之前连接起来。如chr1 TranscriptID条目所示，有两个坐标对(TranscriptID长度)：15和33)，它们应该形成一个序列。基于GeneID的级联是不可能的，因为相同的GeneID (chr11的最后两个条目)的转录is可能不同，这将导致错误的序列。此外，这些序列应该是链特异性的，因此nucleotide.fa.需要反向补充(-)链序列。

我试图获得的核苷酸序列是基于的床上工具getfasta (https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html)。然而，正如在文档中所看到的，所需的输出是不可能的，但是对于我的下游质量控制来说确实是必要的。

示例nucleotide.fa的转换在faTrans中确实很简单。因此，对于如何将nucleotide.fa.转换为origin.bed，我非常感激。您是否知道其他的软件包或转换工具是合适的？

谢谢!

Answer 1

根据您链接的床头工具文档，它似乎会做您想做的事情，使用-s选项反向补充(-)链上的特性，并使用-split选项连接相同文本的外显子。

您首先需要将您的床文件转换为"BED12“格式，每记录一行，当使用-split选项运行时，便携工具getfasta期望这样做。

您可以编写自己的代码来完成此转换，或者使用类似于bed6ToBed12的脚本。下面是一个例子：

awk 'BEGIN{FS=OFS="\t"} {print $1, $2, $3, gensub(/_.*/, "", 1, $4), 0, $5}' origin.bed |\
bed6Tobed12.sh -i stdin > origin.bed12

(前面的AWK调用只是将输入文件转换为正确的BED6格式，其中字符串是第6个字段，并将文本名称提取到第四个字段中。)

在此之后，您可以运行bedtools getfasta来获取核苷酸序列：

bedtools getfasta -fi reference.fa -bed origin.bed12 -s -name -split > nucleotide.fa