寫在前面
要讓從基因組中獲取信息變得像喝水那樣簡單周瞎。 --by zhougengxu
預(yù)先準(zhǔn)備
我以水稻基因組為例，像大家介紹如何提取序列樊销。
基因組文件：水稻基因組整慎，水稻注釋文件，水稻CDS文件(用作檢驗)
工具：samtools围苫，bedtools,awk
1.下載文件
wget https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_genome.fasta.gz
gunzip IRGSP-1.0_genome.fasta.gz
wget https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_2020-12-02.tar.gz
tar xvf IRGSP-1.0_representative_2020-12-02.tar.gz
2.提取基因序列
2.1提取gff文件的所有基因位置,并轉(zhuǎn)換成bed格式
awk '{if($3~/^gene$/)print}' locus.gff > genes.gff
gff2bed <genes.gff> genes.bed
2.2計算染色體長度
samtools faidx IRGSP-1.0_genome.fasta
cut -f 1,2 IRGSP-1.0_genome.fasta.fai > sizes.chr
2.3在基因組序列中抓取基因的序列
bedtools getfasta -s -fi IRGSP-1.0_genome.fasta -bed genes.bed -fo genes.fa -name
打包成腳本如下
#bin/bash
GFF=locus.gff
GENOME=IRGSP-1.0_genome.fasta
awk '{if($3~/^gene$/)print}' $GFF > genes.gff
gff2bed <genes.gff> genes.bed
samtools faidx $GENOME
cut -f 1,2 $GENOME.fai > sizes.chr
bedtools getfasta -s -fi $GENOME -bed genes.bed -fo genes.fa -name
3.提取基因組中的一段序列
比如我只想提取一號染色體上的10000-50000這一段序列裤园。
首先建立索引。
samtools faidx IRGSP-1.0_genome.fasta
然后根據(jù)染色體信息和物理位置直接提取剂府。這里注意拧揽，如何要和目的4連用的話，要修改>后面的值與gff文件一致才可以腺占。
samtools faidx IRGSP-1.0_genome.fasta chr01:10000-50000 > chr1.fa
>chr01:10000-50000
TAAGGAAAGTATAGTCAGGTTTACTTATAATGCAACTTATTGGCATATTAATGCATGTTT
ATTGTTTGCACTCCTGAAGTTGAGTTGGTCAATTTTTTCCATAGATGTTCAATGATAATA
TGCTGAAGATCAATGTAAAAAGATGCGTACGCATGGCTATCAAGCTGAGGAAGAagtata
tttacaaggtaaattgtgataagcatacatgagcaccatccatgtttttcaggtatattt
acaaggtaagttctaaTTTGGTCTTGGGATACAGTTATTGAAGGGGGGCTCGGAATGATG
AAAAAGTGGAAACTGAAGTTGAGAAATCGCCCTGGGGAACTAGATGTTTTTCAACCATGA
ACTTCTCGTAATTATGGGTGATGGAACTTTATATATTCATTAGAAGTCCAAAAAAAATCC
GTATACTCACGGTATGGATTCCCTCTTTCTGTTTATGTTTTGTTAGAGATGGATGACGTA
CATAATCTGCTTGTCTTGTTTCAGGGAATCTGTACATATGCACAGTATACTTGAAAACAA淤袜。。衰伯。铡羡。
4.提取3中對應(yīng)的gff文件
首先提取目的區(qū)域的文件。
less -SN locus.gff | grep 'chr01' | awk '$4>=10000{print$0}' | awk '$5<=50000{print$0}' > trunted.gff
1 chr01 irgsp1_locus gene 11218 12435 . + . ID=Os01g0100200;Name=Os01g0100200;Note=Conserved hypothetica
2 chr01 irgsp1_locus gene 11372 12284 . - . ID=Os01g0100300;Name=Os01g0100300;Note=Cytochrome P450 domai
3 chr01 irgsp1_locus gene 12721 15685 . + . ID=Os01g0100400;Name=Os01g0100400;Note=Similar to Pectineste
4 chr01 irgsp1_locus gene 12808 13978 . - . ID=Os01g0100466;Name=Os01g0100466;Note=Hypothetical protein.
5 chr01 irgsp1_locus gene 16399 20144 . + . ID=Os01g0100500;Name=Os01g0100500;Note=Immunoglobulin-like d
6 chr01 irgsp1_locus gene 22841 26892 . + . ID=Os01g0100600;Name=Os01g0100600;Note=Single-stranded nucle
7 chr01 irgsp1_locus gene 25861 26424 . - . ID=Os01g0100650;Name=Os01g0100650;Note=Hypothetical gene. (O
8 chr01 irgsp1_locus gene 27143 28644 . + . ID=Os01g0100700;Name=Os01g0100700;Note=Similar to 40S riboso
9 chr01 irgsp1_locus gene 29818 34453 . + . ID=Os01g0100800;Name=Os01g0100800;Note=Protein of unknown fu
10 chr01 irgsp1_locus gene 35623 41136 . + . ID=Os01g0100900;Name=Os01g0100900;Note=Sphingosine-1-phospha
其次對目標(biāo)區(qū)域的物理位置進(jìn)行相加減意鲸，便可得出目的區(qū)域中的基因組注釋信息烦周。這里注意，減去的數(shù)值是目的位置的數(shù)值減1怎顾，具體原因就是小學(xué)的切樹理論论矾。這樣才能保證你用截短的注釋文件提取的序列和原文件中提取的序列保持一致。
less -SN trunted.gff | awk '{print$1"\t"$2"\t"$3"\t"$4-9999"\t"$5-9999"\t"$6"\t"$7"\t"$8"\t"$9}' > trunted2.gff
1 chr01 irgsp1_locus gene 1219 2436 . + . ID=Os01g0100200;Name=Os01g0100200;Note=Conserved
2 chr01 irgsp1_locus gene 1373 2285 . - . ID=Os01g0100300;Name=Os01g0100300;Note=Cytochrome
3 chr01 irgsp1_locus gene 2722 5686 . + . ID=Os01g0100400;Name=Os01g0100400;Note=Similar
4 chr01 irgsp1_locus gene 2809 3979 . - . ID=Os01g0100466;Name=Os01g0100466;Note=Hypothetical
5 chr01 irgsp1_locus gene 6400 10145 . + . ID=Os01g0100500;Name=Os01g0100500;Note=Immunoglobulin-like
6 chr01 irgsp1_locus gene 12842 16893 . + . ID=Os01g0100600;Name=Os01g0100600;Note=Single-stranded
7 chr01 irgsp1_locus gene 15862 16425 . - . ID=Os01g0100650;Name=Os01g0100650;Note=Hypothetical
8 chr01 irgsp1_locus gene 17144 18645 . + . ID=Os01g0100700;Name=Os01g0100700;Note=Similar
9 chr01 irgsp1_locus gene 19819 24454 . + . ID=Os01g0100800;Name=Os01g0100800;Note=Protein
10 chr01 irgsp1_locus gene 25624 31137 . + . ID=Os01g0100900;Name=Os01g0100900;Note=Sphingosine-1-phospha
最后得到的trunted2.gff便是我們想要的目標(biāo)3里面基因組文件對應(yīng)的注釋文件了杆勇。
參考鏈接
1.https://www.omicsclass.com/article/1169
2.https://mp.weixin.qq.com/s/0CvJNGM_MJGm-wbOxtczzw
