写在前面
要让从基因组中获取信息变得像喝水那样简单。 --by zhougengxu
预先准备
我以水稻基因组为例，像大家介绍如何提取序列。
基因组文件：水稻基因组，水稻注释文件，水稻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文件的所有基因位置,并转换成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
然后根据染色体信息和物理位置直接提取。这里注意，如何要和目的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中对应的gff文件
首先提取目的区域的文件。
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
其次对目标区域的物理位置进行相加减，便可得出目的区域中的基因组注释信息。这里注意，减去的数值是目的位置的数值减1，具体原因就是小学的切树理论。这样才能保证你用截短的注释文件提取的序列和原文件中提取的序列保持一致。
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便是我们想要的目标3里面基因组文件对应的注释文件了。
参考链接
1.https://www.omicsclass.com/article/1169
2.https://mp.weixin.qq.com/s/0CvJNGM_MJGm-wbOxtczzw
