反补fasta文件中的SOME序列

Question

我读了很多关于反向互补序列的有用文章，但我收到了一个似乎不寻常的请求。我在bash工作，我的stdout中有fasta格式的DNA序列，我想通过管道传递。看似不寻常的一点是，我正试图对其中一些序列进行反向补码，这样输出的所有序列都在同一方向上(用于稍后的多序列比对)。

我的fasta标题以"；C"；或"+&"；。我想对以"；C"；。这里有一个子集：

>chr1:86214203-86220231+
CTGGTGGTACAGCTACATTGTACCATAAAACTTATTCATATTAAAACTTA
TTTATATGTACCTCAAAAGATTAAACTGGGAGATAAGGTGTGGCATTTTT
>chr1:84518073-84524089C
caccttagagataatgaagtatattcagaatgtagaacattctataagac
aactgacccaatatcttttaaaaagtcaatgccatgttaaaaataaaaag

我知道有很多方法可以反向补码，比如：

echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev

和

seqtk seq -r in.fa > out.fa

但我不知道如何只对那些在标头末尾有C的序列执行此操作。我认为awk或sed可能是最好的选择，但我不知道如何真正编码它。我可以用awk获得序列头，比如：

awk '/^>/ { print $0 }'

>chr1:84518073-84524089C
>chr1:86214203-86220231+

但如果有人能帮我弄清楚如何把那个awk语句变成一个"；如果标题中的最后一个字符有C，请执行此操作"那太好了！

已编辑以添加：当我写这篇文章的时候，我太累了，我很抱歉没有包括我想要的输出。以下是我想要输出的内容，使用我的小示例：

>chr1:86214203-86220231+
CTGGTGGTACAGCTACATTGTACCATAAAACTTATTCATATTAAAACTTA
TTTATATGTACCTCAAAAGATTAAACTGGGAGATAAGGTGTGGCATTTTT
>chr1:84518073-84524089C
ctttttatttttaacatggcattgactttttaaaagatattgggtcagtt
gtcttatagaatgttctacattctgaatatacttcattatctctaaggtg

您可以看到以+结尾的序列是不变的，但具有以C结尾的标头的序列是反向补的。

谢谢！

Answer 1

早期的一个答案(由Ed Morton提供)使用自包含的awk过程来选择性地反转注释行后的补码序列，该注释行以"结尾；C"；。尽管我认为这是最好的方法，但我将提供一种可能具有更广泛适用性的替代方法。

这里的过程使用awk的system()函数将从awk中的fasta文件中提取的数据发送到shell，在那里序列可以由现有的用于序列操作的许多shell应用程序中的任何一个进行处理。

我定义了一个awk用户函数来将awk的隔离序列传递给shell。它可以从awk过程的任何部分调用：

function processSeq(s)
{system("echo "" s "" |  tr ACGTacgt TGCAtgca | rev ");}

system函数的自变量是一个字符串，其中包含您要在terminal中键入的命令，以实现所需的结果(在本例中，我使用了问题中提到的一个示例反向补码例程)。需要注意的部分是将出现在shell命令中的引号的正确转义，以及在调用函数时将替换指定给它的序列字符串的变量s。s的值与上面显示的system()的参数中在其前后引用的字符串连接。

隔离所需序列

程序的其余部分介绍如何实现：

"；如果标题中的最后一个字符具有C

在使用shell应用程序之前，awk需要隔离要处理的文件的各个部分。一般而言，awk采用一个或多个pattern/action块，其中只有与给定pattern匹配的records(默认情况下为行)由后续action命令处理。例如，如果pattern/^>/ && /C$/对于该行为真(其中/^>/在一行的开始处寻找">"，而/C$/在同一行的结束处寻找"C")，则以下说明性过程执行打印整行print $0的action

/^>/ && /C$/{ print $0 }

对于当前需要，序列从以>开始、以C结束的任何记录之后的下一个记录(行)开始。引用下一行的一种方法是在遇到C行时设置一个变量(在我的示例中称为line)，并为数值比line变量多一个的记录建立一个稍后的模式。

因为fasta序列可能会延伸到几行，所以我们必须在C标题行之后累积几行连续的行。我通过连接C标题行之后的每一行来实现这一点，直到再次遇到以>开头的记录(或者直到使用END块到达文件末尾)。

为了忽略非C标题行后面的序列行，我使用了一个名为flag的变量，当遇到标题记录时，该变量的值为"do"或"ignore"。

对使用system()命令的自定义函数processSeq()的调用是在C标题action块的开始处进行的，如果变量seq保存累积序列(并且在END块中，对于出现在没有标题行的文件末尾的相关序列)。

测试文件和程序

使用示例fasta的修改版本来测试该过程。它包含一个额外相关的C记录，其中有三条和一条位线，而不是两条，以及一个额外无关的+记录。

seq.fasta:

>chr1:86214203-86220231+
CTGGTGGTACAGCTACATTGTACCATAAAACTTATTCATATTAAAACTTA
TTTATATGTACCTCAAAAGATTAAACTGGGAGATAAGGTGTGGCATTTTT
>chr1:84518073-84524089C
caccttagagataatgaagtatattcagaatgtagaacattctataagac
aactgacccaatatcttttaaaaagtcaatgccatgttaaaaataaaaag
>chr1:86214203-86220231+
CTGGTGGTACAGCTACATTGTACCATAAAACTTATTCATATTAAAACTTA
TTTATATGTACCTCAAAAGATTAAACTGGGAGATAAGGTGTGGCATTTTT
>chranotherC
aatgaagtatattcagaatgtagaacattaactgacccgccatgttaatc
aatatctataagaccttttaaaaagcaccttagagattcaataaagtcag
gaagtatattcagaatgtagaacattaactgactaagaccttttaacatg
gcattgact

程序

awk '
/^>/ && /C$/{
if (length(seq)>0) {processSeq(seq); seq="";} 
line=NR; print $0; flag="do"; next;
} 

/^>/ {line=NR; flag="ignore"}

NR>1 && NR==(line+1) && (flag=="do"){seq=seq $0; line=NR; next} 

function processSeq(s)
{system("echo "" s "" |  tr ACGTacgt TGCAtgca | rev ");}

END { if (length(seq)>0) processSeq(seq);} 
' seq.fasta

输出

>chr1:84518073-84524089C
ctttttatttttaacatggcattgactttttaaaagatattgggtcagttgtcttatagaatgttctacattctgaatatacttcattatctctaaggtg
>chranotherC
agtcaatgccatgttaaaaggtcttagtcagttaatgttctacattctgaatatacttcctgactttattgaatctctaaggtgctttttaaaaggtcttatagatattgattaacatggcgggtcagttaatgttctacattctgaatatacttcatt

在Raspberry Pi 400上使用GNU Awk 5.1.0进行测试。

性能说明

因为调用sytstem()会创建一个子shell，所以这个过程将比自包含的awk过程慢。在现有shell例程可用或难以使用自定义awk例程进行复制的情况下，它可能很有用。

编辑：修改以包括未更改的+记录

这个版本有一些早期块的重复，有一些小的变化，以处理不需要反向补充的行的打印(如果理解了主要解释，这些变化应该是不言自明的)

awk '
/^>/ && /C$/{
if (length(seq)>0 && flag=="do") {processSeq(seq)} else {print seq} seq="";line=NR; print $0; flag="do"; next;
} 

/^>/ {if (length(seq)>0 && flag=="do") {processSeq(seq)} else {print seq} seq=""; print $0; line=NR; flag="ignore"}

NR>1 && NR==(line+1){seq=seq $0; line=NR; next} 
function processSeq(s)
{system("echo "" s "" |  tr ACGTacgt TGCAtgca | rev ");}

END { if (length(seq)>0 && flag=="do") {processSeq(seq)} else {print seq}} 
' seq.fasta

Answer 2

使用任何awk:

$ cat tst.awk
/^>/ {
if ( NR > 1 ) {
prt()
}
head = $0
tail = ""
next
}
{ tail = ( tail == "" ? "" : tail ORS ) $0 }
END { prt() }
function prt(    type) {
type = substr(head,length(head),1)
tail = ( type == "C" ? rev( tr( tail, "ACGTacgt TGCAtgca" ) ) : tail )
print head ORS tail
}
function tr(oldStr,trStr,   i,lgth,char,newStr) {
if ( !_trSeen[trStr]++ ) {
lgth = (length(trStr) - 1) / 2
for ( i=1; i<=lgth; i++ ) {
_trMap[trStr,substr(trStr,i,1)] = substr(trStr,lgth+1+i,1)
}
}
lgth = length(oldStr)
for (i=1; i<=lgth; i++) {
char = substr(oldStr,i,1)
newStr = newStr ( (trStr,char) in _trMap ? _trMap[trStr,char] : char )
}
return newStr
}
function rev(oldStr,    i,lgth,char,newStr) {
lgth = length(oldStr)
for ( i=1; i<=lgth; i++ ) {
char = substr(oldStr,i,1)
newStr = char newStr
}
return newStr
}

$ awk -f tst.awk file
>chr1:86214203-86220231+
CTGGTGGTACAGCTACATTGTACCATAAAACTTATTCATATTAAAACTTA
TTTATATGTACCTCAAAAGATTAAACTGGGAGATAAGGTGTGGCATTTTT
>chr1:84518073-84524089C
ctttttatttttaacatggcattgactttttaaaagatattgggtcagtt
gtcttatagaatgttctacattctgaatatacttcattatctctaaggtg

Answer 3

这可能对你有用(GNU sed)：

sed -nE ':a;p;/^>.*C$/!b
:b;n;/^>/ba;s/^/n/;y/ACGTacgt/TGCAtgca/
:c;tc;/n$/{s///p;bb};s/(.*)n(.)/21n/;tc' file

打印当前行，然后检查。

如果线路不是以>开始，而以C结束，则拉出并重复。

否则，获取下一行，如果它以>开头，则重复上一行。

否则，插入一条换行符(反转该行时用作轴心点)，使用平移命令补充该行的代码。然后开始逐字符反转行，直到插入的换行符到达行的末尾。

删除换行符，打印结果并重复上面的行。

注意：n命令将在读取最后一行之后执行脚本时终止脚本。

---

由于OP修改了输出，另一种解决方案是当整个序列被补码然后反转时。以下是另一个我认为符合这些标准的解决方案。

sed -nE ':a;p;/^>.*C$/!b
:b;n;/^>/!{H;$!bb};x;y/ACGTacgtn/TGCAtgca%/;s/%/n/
:c;tc;s/n$//;td;s/(.*)n(.)/21n/;tc
:d;y/%/n/;p;z;x;$!ba' file