提取序列是生物資訊分析中常見的一個操作，也是學習生物資訊程式設計的入門操作。通常是給定基因ID，然後從一個大的資料集裡面提取出匹配ID的序列，包含匹配的序列ID和序列資訊，類似於Excel中的Vlookup，但是這裡需要一個包含序列ID的列表以及一個包含序列的fasta格式檔案。如果不會程式設計該如何提取呢，今天我們就介紹一些方法。
例如這裡有五條序列，我們需要根據基因ID，提取出gene3和gene5的內容。
>
gene1
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTCTCTGACAGCAGC
TTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAA
TATAGGCATAGCGCACAGACAGATAAAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACC
>
gene2
ATTACCACCACCATCACCACCACCATCACCATTACCATTACCACAGGTAACGGTGCGGGCTGACGCGTAC
AGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTCGACCAAAGGTAACGAGGTAACA
>
gene3
ACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGGGTTGCCGATA
TTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCA
CCTGGTGGCGATGATTGAAAAAACCATTAGCGGCCAGGATGCTTTACCCAATATCAGCGATGCCGAACGT
ATTTTTGCCGAACTTCTGACGGGACTCGCCGCCGCCCAGCCGGGATTCCCGCTGGCGCAATTGAAAACTT
>
gene4
TCGTCGACCAGGAATTTGCCCAAATAAAACATGTCCTGCATGGCATTAGTTTGTTAGGGCAGTGCCCGGA
TAGCATTAACGCTGCGCTGATTTGCCGTGGCGAGAAAATGTCGATCGCCATTATGGCCGGCGTATTAGAA
>
gene5
GCGCGCGGTCACAACGTTACCGTTATCGATCCGGTCGAAAAACTGCTGGCAGTGGGGCATTACCTCGAAT
CTACTGTCGATATTGCAGAGTCCACCCGCCGTATTGCGGCAAGTCGTATTCCGGCTGATCACATGGTGCT
GATGGCAGGTTTCACCGCCGGTAATGAAAAAGGCGAACTGGTGGTACTTGGACGCAACGGTTCCGACTAC
TCCGCGGCGGTGCTGGCTGCCTGTTTACGCGCCGATTGTTGCGAGATTTGGACGGACGTTGACGGGGTAT
原始方法
直接用文字編輯器開啟，然後直接Ctrl+C複製，Ctrl+V貼上，不過這樣處理稍微大一點的資料，不僅效率低下，而且很容易出錯，直至就會把人搞奔潰，我們堅決抵制這種做法。
sed
利用sed提取，我們首先使用less -N檢視一下這兩條序列對應的行號，然後就可以利用sed輸出任意行的內容了。
less
-N gene。fna
sed -n
‘8，12p；16，20p’
gene。fna
awk
awk也可以輸出固定行的內容，或者匹配固定行的內容。其中NR代表行號，number row。
awk
‘NR>=8 && NR<12 {print}’
gene。fna
awk
‘NR>=16 && NR<20 {print}’
gene。fna
grep
grep可以用於匹配ID，-A選項可以設定輸出匹配後的幾行，我們首先利用seqtk工具將fasta序列都格式化為兩行一個單位。
seqtk
seq
-l
0
gene
。fna
>
gene
。fa
grep
-A
1 “
gene
［3|5］
”
gene
。fa
samtools
以上幾種方法處理一些小序列還行，如果一次有100個基因ID，也不太方便。
這種情況下強烈推薦samtools工具。利用samtools建立索引，就可以完成快速提取序列。
#首先為利用faidx為fasta檔案建立索引
samtools
faidx gene。fna
#建立索引之後就可以快速提取了
samtools faidx gene。fna gene3 gene5
程式設計
其實用程式設計也很容易實現，無論是perl還是python都不難，首先將ID和序列儲存為一個雜湊或者字典型的資料結構，然後就可以很方便的利用ID進行查找了。
#將基因
ID
寫入到
gene
。list
檔案中，每行一個基因
ID
perl
get_seq_bylist
。pl
gene
。list
gene
。fna
