genblastg和genewise鉴定缺失基因
2023-11-07 本文已影响0人
多啦A梦的时光机_648d
基于组装好的基因组，例如某一个基因在该基因组中无注释，需检查是否为注释错误，还是该基因本身缺失。
软件：genblastg和genewise
genblastg下载地址：https://github.com/thecgs/genblastG_extension
1.下载目标基因的蛋白序列，以拟南芥CYP711A1基因为例，到官网下载https://www.arabidopsis.org/servlets/TairObject?id=34142&type=locus，保存为CYP711A1.faa
2.基因组序列文件：final_assembly_descending.fa
3. 利用genblastg在基因组序列中查找，找到基因组中CYP711A1基因的候选区域。并整理成bed文件
export PATH=/home/lx_sky6/yt/soft/miniconda3/envs/genblastg/bin/:$PATH
python3 /home/lx_sky6/yt/soft/genblastG_extension-master/run_formatdb_nucl.py final_assembly_descending.fa
python3 /home/lx_sky6/yt/soft/genblastG_extension-master/run_genblastG.py -g final_assembly_descending.fa -i CYP711A1.faa -o out.gff3
###注意，脚本中-c默认为0.8，可自行修改阈值
```{tmp.bed}
Chr7 44465688 44468002
Chr7 44490184 44504841
Chr7 54220158 54222085
Chr7 54225358 54227310
4. 向该区域两端延伸500bp,并利用bedtools提取序列
cat tmp.bed |awk '{print $1"\t"$2-500"\t"$3+500}' >ex500.bed
Chr7 44465188 44468502
Chr7 44489684 44505341
Chr7 54219658 54222585
Chr7 54224858 54227810
bedtools getfasta -fi final_assembly_descending.fa -bed 500.bed -fo ex500.fa
5. 然后利用拟南芥CYP711A1.faa为参考序列，利用genewise软件简单CYP711A1基因结构。由于genewise每次只能做一条序列，所以将ex500.fa序列分成4份。
seqkit split ex500.fa -p 4
cd ~/yt/test/genewise/ex500.fa.split
genewise ../CYP711A1.faa ex500.part_001.fa -both -gff >1.out
genewise ../CYP711A1.faa ex500.part_002.fa -both -gff >2.out
genewise ../CYP711A1.faa ex500.part_003.fa -both -gff >3.out
genewise ../CYP711A1.faa ex500.part_004.fa -both -gff >4.out
6. 剔除内含子，提取蛋白编码序列
Chr7:44465188-44468502 GeneWise gene 503 2595 544.59 + . Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise cds 503 750 0.00 + 0 Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise intron 751 911 0.00 + . Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise cds 912 1064 0.00 + 1 Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise intron 1065 1213 0.00 + . Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise gene 503 2595 544.59 + . ID=CYP711A1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise mRNA 503 2595 544.59 + . ID=CYP711A1.1;Parent=CYP711A1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise exon 503 750 0.00 + 0 Parent=CYP711A1.1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise cds 503 750 0.00 + 0 Parent=CYP711A1.1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise exon 912 1064 0.00 + 1 Parent=CYP711A1.1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise cds 912 1064 0.00 + 1 Parent=CYP711A1.1;Chr7:44465188-44468502-genewise-prediction-1
Chr7:44465188-44468502 GeneWise exon 1214 1997 0.00 + 1 Parent=CYP711A1.1;Chr7:44465188-44468502-genewise-prediction-1
然后利用gffread提取cds和pep
gffread 1.gff -g ex500.part_001.fa -x cds.fa -p pep.fa
>CYP711A1.1
NETLVRSMAALDWDALTTKSASTAFTMLAILAGVLVYFYAPYWGVRRVPGPPALPILGHLHLIVKYGTSI
FSVLADQYGPVFRFHMGRQPLVIIADPELCREVGIKKFKDFSDRSMPYSISTIAVHANGYPFSREPRRSA
TRSILLSIYQSSFVANLLPTMQAYIESAAENLNIKDEESIDFFTLSHKMTTDIIGKAAFGVNFGLSRCLN
SAALDTSDIQIISEKDASEFIDKHTYATTCLKMDLAGSFSMIVGQILPLFQKPLYHFLRRIPGTMDYKIF
TSHEDLGDRINKIVDSRIQHGDTGSKDFLSLILKATESETENKRFFTRDYVGALAYEHIIGGSTSTAVTL
TCIIYLVAGHPEVQKKLLNEIDAFGPIDKVPTPEDLQNNFPYLDLVVKEGMRYYTLSPITAREAGRDVHI
GGYFIPKGTWVWPALGVLAKDPKHFPDPEAFKPERFDPNGEEEKNRHPYAYIPFGIGPRACLGQKFALQV
LKLSLIHQYRRYVFARDGPIEFQYGVFLTFKTEVKFKVTKRC
和拟南芥序列进行比对
CLUSTAL format alignment by MAFFT L-INS-i (v7.520)
AT2G26170.1 MKTQHQWWEVLD-PFLTQHEALIAFLTFAAVV-IVIYLYRPSWSVCNVPGPTAMPLVGHL
AT2G26170.2 ------------------------------------------------------------
CYP711A1.1 NETLVRSMAALDWDALTTKSASTAFTMLAILAGVLVYFYAPYWGVRRVPGPPALPILGHL
AT2G26170.1 PLMAKYGPDVFSVLAKQYGPIFRFQMGRQPLIIIAEAELCREVGIKKFKDLPNRSIPSPI
AT2G26170.2 -------------------------MGRQPLIIIAEAELCREVGIKKFKDLPNRSIPSPI
CYP711A1.1 HLIVKYGTSIFSVLADQYGPVFRFHMGRQPLVIIADPELCREVGIKKFKDFSDRSMPYSI
******:***:.*************:.:**:* .*
AT2G26170.1 SASPLHKKGLFFTRDKRWSKMRNTILSLYQPSHLTSLIPTMHSFITSATHNLDSKPRD-I
AT2G26170.2 SASPLHKKGLFFTRDKRWSKMRNTILSLYQPSHLTSLIPTMHSFITSATHNLDSKPRD-I
CYP711A1.1 STIAVHANGYPFSREPRRSATRSILLSIYQSSFVANLLPTMQAYIESAAENLNIKDEESI
*: .:* :* *:*: * * *. :**:**.*.::.*:***:::* **:.**: * .: *
AT2G26170.1 VFSNLFLKLTTDIIGQAAFGVDFGL-----------SGKKPIKDVEVTDFINQHVYSTTQ
AT2G26170.2 VFSNLFLKLTTDIIGQAAFGVDFGL-----------SGKKPIKDVEVTDFINQHVYSTTQ
CYP711A1.1 DFFTLSHKMTTDIIGKAAFGVNFGLSRCLNSAALDTSDIQIISEKDASEFIDKHTYATTC
* .* *:******:*****:*** *. : *.: :.::**::*.*:**
AT2G26170.1 LKMDLSGSLSIILGLLIPILQEPFRQVLKRIPGTMDWRVEKTNARLSGQLNEIVSKRAKE
AT2G26170.2 LKMDLSGSLSIILGLLIPILQEPFRQVLKRIPGTMDWRVEKTNARLSGQLNEIVSKRAKE
CYP711A1.1 LKMDLAGSFSMIVGQILPLFQKPLYHFLRRIPGTMDYKIFTSHEDLGDRINKIVDSRIQH
CYP711A1.1 基因其实是存在于基因组中，注释漏掉了，因此需要手动添加注释到gff中.
修改为chr上正确的位置
44466106 44466252
44466402 44467185
44467277 44467372
44467469 44467783
Chr7 EVM gene 44466106 44467783 544.59 + . ID=evm.TU.Chr7.3604;Name=EVM%20prediction%20Chr7.3604
Chr7 EVM mRNA 44466106 44467783 544.59 + . ID=evm.model.Chr7.3604;Parent=evm.TU.Chr7.3604;Name=EVM%20prediction%20Chr7.3604
Chr7 EVM exon 44466106 44466252 0.00 + 1 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM cds 44466106 44466252 0.00 + 1 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM exon 44466402 44467185 0.00 + 1 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM cds 44466402 44467185 0.00 + 1 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM exon 44467277 44467372 0.00 + 0 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM cds 44467277 44467372 0.00 + 0 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM exon 44467469 44467783 0.00 + 0 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
Chr7 EVM cds 44467469 44467783 0.00 + 0 Parent=evm.model.Chr7.3604;Parent=evm.model.Chr7.3604
