我读了很多关于反向互补序列的有用文章，但我收到了一个似乎不寻常的请求。我在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结尾的标头的序列是反向补的。
谢谢！
早期的一个答案(由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
使用任何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
这可能对你有用(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