目標(biāo)：下載Synechococcus elongatus UTEX 2973 的基因組注釋文件，統(tǒng)計其中染色體序列(CP006471.1)前10kb有幾個基因？
1.搜索基因
NCBI的Assembly數(shù)據(jù)庫是一個綜合基因組數(shù)據(jù)庫，首先打開NCBI，將搜索框勾選為Assembly，輸入需要搜索的基因名稱，搜索目的基因組，搜索結(jié)果如圖1.1，在頁面左邊可以看到它的基本信息，包括Organism name ，accession number等，在頁面右邊可以看到數(shù)據(jù)庫來源。
2.在瀏覽器中復(fù)制文件路徑
GenBank和Refseq數(shù)據(jù)庫的區(qū)別在于Refseq源于Genbank，提供非冗余序列信息，這里直接點擊1.1右側(cè)FTP directory for GenBank assembly就可以，點擊后可以看到需要的各種數(shù)據(jù)文件的鏈接，如圖2.1，此次下載的是基因組的注釋文件，也就是以gff為后綴的文件。
3.使用wget命令下載并解壓
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/817/325/GCA_000817325.1_ASM81732v1/GCA_000817325.1_ASM81732v1_genomic.gff.gz
接下來使用gunzip解壓gz文件
gunzip GCA_000817325.1_ASM81732v1_genomic.gff.gz
less GCA_000817325.1_ASM81732v1_genomic.gff
4.grep抓取所需要的信息
對于gff文件，第一列通常代表染色體序列信息，第二列代表基因結(jié)構(gòu)來源，第三列代表的是類型即gene、mRNA或其他，第四、五列分別是對應(yīng)片段的起始和終止位置。
統(tǒng)計特定染色體序列(CP006471.1)前10kb的基因數(shù)目可以用grep命令抓取,也就是該基因的終止位置(第五列)小于10k。
grep '^CP006471.1' GCA_000817325.1_ASM81732v1_genomic.gff |awk -v FS="\t" -v OFS="\t" '{if($5<10000) {print $5}}'|sort|uniq|wc -l
最終可以得到染色體序列(CP006471.1)前10kb共有18個基因
關(guān)于上述代碼的具體解釋：
1.grep 用于查找文件里符合條件的字符串，^用于模式的最左側(cè)，這里'^CP006471.1'即匹配以CP006471.1開頭的內(nèi)容。
grep完整的語法結(jié)構(gòu)
grep [option] [pattern]
2.awk傾向于一行當(dāng)中分成數(shù)個字段來處理,其中參數(shù)-v表示定義或修改一個awk內(nèi)部變量，awk內(nèi)置變量FS表示輸入字段分隔符，默認(rèn)為空格符；內(nèi)置變量OFS表示輸出字段分隔符，默認(rèn)也為空格符。
此外awk中的-F參數(shù)也可以用來指定分隔符。
awk語法結(jié)構(gòu)
awk [option] 'pattern {action}'
3.sort|uniq|wc -l用以統(tǒng)計
sort以行為單位對文件進行排序
uniq刪除相鄰的重復(fù)的行并寫到標(biāo)準(zhǔn)輸出
wc統(tǒng)計指定文本文件的行數(shù)、字?jǐn)?shù)、字符數(shù)，其中-l參數(shù)表示統(tǒng)計行數(shù)
經(jīng)過上述操作我們就完成了對基因組注釋文件的下載與基因數(shù)量的統(tǒng)計。
