写在前面
看到一个问题，如何筛选fasta文件中基因序列长度，当时脑子里想到了最近正在学的东西。如何在linux下一行指令筛选fasta文件中序列长度低于200的序列。
cat read_1.fa | paste - - | awk 'length($2) <=200 {print}' | tr "\t" "\n" > low_seq.fasta
命令都是连着的 ：cat用于读取文件 paste 是把前两行放在一起用制表符分隔，awk判断并输出长度低于200的序列最后用tr把空格改成换行符 输出到一个新的文件 问题似乎解决了。 但真的这么容易么？实际操作的过程却发现不太好使。
文件是gzip格式的，cat命令没有办法正确识别啊……那就传入解压命令
cat <(gunzip -c cdna.fa.gz) | paste - - | awk 'length($2) <=200 {print}' | tr "\t" "\n" | less
#结果
>TraesCS5D02G061900.1 cdna chromosome:IWGSC:5D:58143379:58149665:1 gene:TraesCS5D02G061900 gene_biotype:nontranslating_CDS transcript_biotype:nontranslating_CDS
ATGCAGCACTTCTTCCTCTTCCTCTTGCTCCTCGTTTTCTTGCTTCCCCATGAAACCAGC
TACTCCGCGGCCGCAGCTTCCGGCGGCGGTGACAGGTGCAGCCGCCGGTGTGGCGGTGCC
ACCGTCCCATACCCTCTCGGATTCTCCCCCGGCTGCCCCATCGCCCTGTGGTGCGACGCC
AGCATATCCACGCCGATCCTCCCATACATAGGAGAGAACGGCACCACGTACCGCGTGATA
GCCTTCAACTCCACCATCTCCACCGTCGTTGTCGGCCTCCCGCCGTCGTGCAGCCGCAGC
GTCCCCGAGGCCAGGAGGGTGCTCTCCGGGGGCAACTACGGCGTGTCGTCTCGCACCGGT
CTTTTCCTCCGCGGCGGCTGCGGCGGGGTCAACGCGTCGGCGGGCGCTGGATGCAGCGTG
CCCGTGGGCGTCATGTCCAGGCTGCTCCGCACGGCGCAGTGCGTCGGCATCAGCAACGAC
ACCTCGCCCGCCGCTGTGGCGTGCGTCGCTTCCGCTGCGCCAAACTCCACCGCGGCGGTG
CATGGCCAGTTTCTGCGGTGGGAGAAGGTCGAGAAGCAAAAGTGCGACGACGTGCTGACC
TCCACGGTGTTCGTGCTCACGGCGGAGGGGACGGCTACGCTGGAGTTCGGCGTGGCCGAG
CTCGGCTGGTGGCTCAACGGGACGTGCGGCGCCGGCGGGGGCGAGCGTTGCGCGGCGAAC
GCGTCGTGCACCGACGTAAAAACGCCTAGCGGGGAGTTGGGGCACCGGTGCGGGTGCGTG
GCGGGGATGGAGGGCGACGGCTTCTTGGCCGGCGACGGATGCTACCTCGGTGGATCCCCT
CATGATGTTGTCTCAAGCTCCAAGGTGCGGCAAATATTCATAGTTATTGGCTTGATATCT
TCTCTTGGCATCACTGTCTTCTTCTTCATCTATGGCCGGAAAAGGCGCAATGCCCTGATG
AAGACAACAAAACAACTACCCAAAAAGGCAGCCATATTCTTCAATGGGGAAGTCATGGAG
GACGAGCTAGAGCAAGGAGCCAGCGGCCCCCGGCGATTTTCCTACAGTGAGCTCGCCGTC
GCCACTGATAACTTTTCAAATGACATGGCACTCGGGAGAGGAGGCTTCGGGTATGTATAT
CAAGGGTTTATTAGTGACATGGACCGTGAAGTGGCTGTCAAGAGGGTGTCTGAGACCTCT
CGACAGGGTTGGAAGGAGTTCGTCTCTGAGGTGCGGATCATTAGTCGGCTCAGGCATCGA
AACCTTGTGCAGCTCATAGGCTGGTGCCACGGTGGAGACGAGCTTCTCCTCGTCTACGAC
CTGATGCACAATGGCAGCCTCGACACTCATTTGTATAGATCCGACCATGGGCTGACATGG
CTGATAAGAGGCGGAGCAACGTGTGGTGCACAGGGACATCAAGCCTAGCAACATCATGTT
AGACGCCTCCTTCACAGCCAAGCTCGGTGACTTCGGGCTCGCGAGGCTCATCAACGACGG
ACGGAGATCGCACACGACCGGCATAGCTGGCACGATGGGGTACATCGATCCAGAGAGCAT
GCTGGCCGGCAGGGCGAGCATCGAGTCAGACGTGTATAGCTTCGGGGTCGTCCTTCTCGA
GGTTGCATGTGGGCGGCGACCAGCTTTGGTCCAGGAAGACGGTGATGTCGTCCACCTGGT
CCAATGGGTGTGGGACTTGTATGGCAGAGGATCAATCTTTGGCGCTGTCGACAAGCGACT
TAGGGGGGAGTTTGACGGCACGGAGATGGAGCGTGTTATGGTCGCGGGGCTCTGGTGCGC
GCACCCTGACCGTGGCATGCGACCATCCATCAGGCAAGTGGTGAACATGCTGCGGTCCGA
AGCGCCATTACCAAGTCTCCCGGCGAGGATGCCGGTGGCGACCTATGGGCCGCCGACTAA
TCATTCGAGTTTTGGAACTTTGCTGATGACCAGTGTCAATGGCCGGTGA
>TraesCS7A02G035884.1 cdna chromosome:IWGSC:7A:15957212:15957450:-1 gene:TraesCS7A02G035884 gene_biotype:nontranslating_CDS transcript_biotype:nontranslating_CDS
ATGAAGACCATGTTCATCCTCGCTCTCCTCACCTTCACGGCGACCAGCGCAGTTGCGCAG
CTGGACACTACCTGTAGCCAGGGCTACGGGCAATGCCAGCAGCAGTCGCAGCCTTAGATG
AACACATGCGCTGCTTTCCTGCAGCAGTGAGCCAGACACCATATGTCCAGTCACAGATGT
GGCAGGCAAGCAGTTGCCAGTTGATGCGGCAACAGTGCTGCCAGCTGACACGTCTTCAT
这跟说好的不一样啊！第一个明显超过200了！
仔细看了一眼原数据
首先在名字那一栏多了好多基因相关的信息，例如基因位置，基因名称，基因类型，转路本类型……其次序列中间也有换行符……所以首先需要标准化一下fasta文件格式。
去掉换行符的命令为：
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n } ' <(gunzip -c cdna.fa.gz) |less
#结果
>TraesCS7A02G035884.1 cdna chromosome:IWGSC:7A:15957212:15957450:-1 gene:TraesCS7A02G035884 gene_biotype:nontranslating_CDS transcript_biotype:nontranslating_CDS
ATGAAGACCATGTTCATCCTCGCTCTCCTCACCTTCACGGCGACCAGCGCAGTTGCGCAGCTGGACACTACCTGTAGCCAGGGCTACGGGCAATGCCAGCAGCAGTCGCAGCCTTAGATGAACACATGCGCTGCTTTCCTGCAGCAGTGAGCCAGACACCATATGTCCAGTCACAGATGTGGCAGGCAAGCAGTTGCCAGTTGATGCGGCAACAGTGCTGCCAGCTGACACGTCTTCAT
>TraesCS1D02G008600.1 cdna chromosome:IWGSC:1D:4341746:4342864:1 gene:TraesCS1D02G008600 gene_biotype:nontranslating_CDS transcript_biotype:nontranslating_CDS
ATGAAGACCTTCCTCATCTTTGCCCTCCTCGCCGTTGCGGCGACAAGTGCCATTGCACAAATGGAGACTAGCCACATCCCTGGCTTGGAGAAACCATCGCAACAACAACCATTACCACTACAACAAATATTATGGTACCACCAACAGCAACCCATCCAACAACAACCACAACCATTTCCACAACAGCCACCATGTTCACAGCAACAACAACCACCATTATCGCAGCAACAACAACCACCATTTTCACAGTAACAACCACCATTCTCGCAGCAAGAACTACCAGTTCTACCGCAACAACCACCATTTTCGCAGCAACAACAACCATTCCCGCAGCAACAACAACCACTTCTACCGCAACAACCCCCATTTTCACAACAACGACCACCATTTTCTCAGCAGCAGCAACAACCAGTTCTACCGCAACAACCACCATTTACGCAACAACAGCAACAACAACAACAACAACAACAACAACAACAACAACAACAACAACAACAACCAATTCTACCGCAACAACCACCTTTTTCGCAACACCAACAACCGGTTCTGCCGCAACAACAAATACCATATGTTCAGCCATCTATCTTGCAGCAGCTAAACCCATGCAAGGTATTCCTCCAGCAGCAATGCAGCCCTGTGGCAATGCCACAAAGTCTTGCTAGGTCACAAATGTTGTGGCAGAGCAGTTGCCATGTGATGCAGCAACAATGTTGCCAGCAGCTGCCGCGAATCCCCGAACAATCACGCTATGATGCAATCCGTGCCATCATCTACTCGATCGTCCTACAAGAACAACAACATGGTCAGGGTTTCAACCAACCTCATCAGCAACAACCCCAACAGTCGGTCCAAGGTGTCTCCCAACCCCAACAACAACAGAAGCAGCTCGGACAGTGTTTTTTCCAACGACCTCAACAACAACAACTGGGTCAATGGCCTCAACAACAACAGGTACCACAGGGTACCTTGTTGCAGCCACACCAAATAGCTCAACTTGAGTTGATGACTTCCATTGCACTCCGTACCCTGCCAATGATGTGCAGTGTCAACGTGCCGGTGTACGGCACCACCACTAGTGTGCCATTCGGCGTTGGCACCCAAGTTGGTGCCTACTGATAA
总的命令就是不是以>开头的就去掉"\n"，是的话就原样输出，最后全部输出。
之后再把>后的基因名做一个简化，不做的话后续步骤中会出问题因为信息不一定等长……当然你也可以选择把第一行里的空格全部换成"" 只需要加入命令 sed 's/ //g'。不展开说了
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n } ' <(gunzip -c cdna.fa.gz) | awk '{if($0~/>/){sub(/\..*/, "", $0); print$0}else{print$0}}' | less
#结果
>TraesCS5D02G061900
ATGCAGCACTTCTTCCTCTTCCTCTTGCTCCTCGTTTTCTTGCTTCCCCATGAAACCAGCTACTCCGCGGCCGCAGCTTCCGGCGGCGGTGACAGGTGCAGCCGCCGGTGTGGCGGTGCCACCGTCCCATACCCTCTCGGATTCTCCCCCGGCTGCCCCATCGCCCTGTGGTGCGACGCCAGCATATCCACGCCGATCCTCCCATACATAGGAGAGAACGGCACCACGTACCGCGTGATAGCCTTCAACTCCACCATCTCCACCGTCGTTGTCGGCCTCCCGCCGTCGTGCAGCCGCAGCGTCCCCGAGGCCAGGAGGGTGCTCTCCGGGGGCAACTACGGCGTGTCGTCTCGCACCGGTCTTTTCCTCCGCGGCGGCTGCGGCGGGGTCAACGCGTCGGCGGGCGCTGGATGCAGCGTGCCCGTGGGCGTCATGTCCAGGCTGCTCCGCACGGCGCAGTGCGTCGGCATCAGCAACGACACCTCGCCCGCCGCTGTGGCGTGCGTCGCTTCCGCTGCGCCAAACTCCACCGCGGCGGTGCATGGCCAGTTTCTGCGGTGGGAGAAGGTCGAGAAGCAAAAGTGCGACGACGTGCTGACCTCCACGGTGTTCGTGCTCACGGCGGAGGGGACGGCTACGCTGGAGTTCGGCGTGGCCGAGCTCGGCTGGTGGCTCAACGGGACGTGCGGCGCCGGCGGGGGCGAGCGTTGCGCGGCGAACGCGTCGTGCACCGACGTAAAAACGCCTAGCGGGGAGTTGGGGCACCGGTGCGGGTGCGTGGCGGGGATGGAGGGCGACGGCTTCTTGGCCGGCGACGGATGCTACCTCGGTGGATCCCCTCATGATGTTGTCTCAAGCTCCAAGGTGCGGCAAATATTCATAGTTATTGGCTTGATATCTTCTCTTGGCATCACTGTCTTCTTCTTCATCTATGGCCGGAAAAGGCGCAATGCCCTGATGAAGACAACAAAACAACTACCCAAAAAGGCAGCCATATTCTTCAATGGGGAAGTCATGGAGGACGAGCTAGAGCAAGGAGCCAGCGGCCCCCGGCGATTTTCCTACAGTGAGCTCGCCGTCGCCACTGATAACTTTTCAAATGACATGGCACTCGGGAGAGGAGGCTTCGGGTATGTATATCAAGGGTTTATTAGTGACATGGACCGTGAAGTGGCTGTCAAGAGGGTGTCTGAGACCTCTCGACAGGGTTGGAAGGAGTTCGTCTCTGAGGTGCGGATCATTAGTCGGCTCAGGCATCGAAACCTTGTGCAGCTCATAGGCTGGTGCCACGGTGGAGACGAGCTTCTCCTCGTCTACGACCTGATGCACAATGGCAGCCTCGACACTCATTTGTATAGATCCGACCATGGGCTGACATGGCTGATAAGAGGCGGAGCAACGTGTGGTGCACAGGGACATCAAGCCTAGCAACATCATGTTAGACGCCTCCTTCACAGCCAAGCTCGGTGACTTCGGGCTCGCGAGGCTCATCAACGACGGACGGAGATCGCACACGACCGGCATAGCTGGCACGATGGGGTACATCGATCCAGAGAGCATGCTGGCCGGCAGGGCGAGCATCGAGTCAGACGTGTATAGCTTCGGGGTCGTCCTTCTCGAGGTTGCATGTGGGCGGCGACCAGCTTTGGTCCAGGAAGACGGTGATGTCGTCCACCTGGTCCAATGGGTGTGGGACTTGTATGGCAGAGGATCAATCTTTGGCGCTGTCGACAAGCGACTTAGGGGGGAGTTTGACGGCACGGAGATGGAGCGTGTTATGGTCGCGGGGCTCTGGTGCGCGCACCCTGACCGTGGCATGCGACCATCCATCAGGCAAGTGGTGAACATGCTGCGGTCCGAAGCGCCATTACCAAGTCTCCCGGCGAGGATGCCGGTGGCGACCTATGGGCCGCCGACTAATCATTCGAGTTTTGGAACTTTGCTGATGACCAGTGTCAATGGCCGGTGA
好！很有精神！
最后再套上原来的操作命令
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n } ' <(gunzip -c cdna.fa.gz) | awk '{if($0~/>/){sub(/\..*/, "", $0); print$0}else{print$0}}' | paste - - | awk 'length($2) <=200 {print}' | tr "\t" "\n"| less
#结果
>TraesCS4A02G403700
ATGGTCAAGCATTTACCTAAGTTATGGAAATTTCATAATTCATCAGTTATTTTAGATAATACTAAAAAACAAGATAGAGTAAGATACAGAAGTGTAAAAGAAAGTGGATTCATCAGATGA
>TraesCS1A02G233300
ATGGAGAACGGACAGCCGCAGCCTCCTCCAGGCTACCCGACGGTGGACCCGAAACAGCAGGCAGGCGGCAGGAGGAGATGCTGCTGCGGGCCGTGGCGCCGCAGGACCAAGCAGCGTGGGCGTGGGGAGACCAGCTTCATCGAAGGATGCCTTGCC
GCGCTCTGTTGCTGCTGGCTTTGCGAGCTGTGCTGCGATTGA
>TraesCS5D02G371900
ATGATTCCTCTGGAATTACGACCTTTACCACAACGGTGTCGTCCATGGATGAAATTATTTCGTGGATTGGATTTCACTTGCCTGTCTACGGTTCCCTTGCGTGTGCTCGGGATAGGTGTTTTGTATAAATGTTTCGCCGTATTATTAAGTATTCTC
CTTTAG
#如果需要保存
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n } ' <(gunzip -c cdna.fa.gz) | awk '{if($0~/>/){sub(/\..*/, "", $0); print$0}else{print$0}}' | paste - - | awk 'length($2) <=200 {print}' | tr "\t" "\n" > output.fa
问题似乎得到了解决，孩子欣慰的笑了。回想整个过程其实挺坎坷的，首先你要查看文件，解压缩，去掉多余的回车，简化ID，判断序列长度最后输出。这是我目前会的办法，我相信老手可以用更简单的一行代码搞定所有问题。但我嘛，还在路上，边走边学习吧。
其实可以更简单的
其实如果用TBtools提取可以完全不管fasta格式也不需要考虑命令，只需要四步走：
1.提取所有序列的长度
长度
2.提取长度小于200的序列ID
ID
-
注意这里的ID后面还是有序列长度数值的，你可以拖入excel（注意是讲文件拖入）获得两行，然后复制ID即可。为了操作连贯性这里直接在TBtools中操作。
1
文件直接拖入打开，选择第一行任意位置，然后selected done
2
设置好输出目录点确定即可
3
-
序列提取
最后，这套操作的应用性其实非常强。熟悉这几个操作你能干出很多事。祝 磕盐顺利……
网友评论