05鉴定duplicate gene
blastp比对
首先使用整个基因组的蛋白序列进行比对
-evalue 1e-10
筛选指标-outfmt 6
下一步 MCScanX 要使用的格式-num_threads 20
线程数目
如果存在多个转录本的情况下默认用第一个转录本代表这个基因
提取第一个转录本代表这个基因的氨基酸序列
awk 'NR==1&&$1~/^>/{print $0}NR>=2&&$1~/^>/{print "\n"$0}$1~/^[^>]/{printf $0}' Ghirsutum_gene_peptide.fasta|grep -E "Ghir_A[^\.]*\.1$" -A 1|sed -e 's/-//g' -e 's/\.1//g'
## 建库
makeblastdb -dbtype 'prot' -parse_seqids -in Dt_peptide.fa -out blastpDB/Dt
# 比对
blastp -query At_peptide.fa -db At -evalue 1e-10 -num_threads 20 -outfmt 6 -out At_vs_At.blast
## 提交LSF任务
bsub -J Dt_Dt -q "normal" -n 20 -R span[hosts=1] -e Dt_Dt.err -o Dt_Dt.out "blastp -query ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Dt_peptide.fa -db ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/blastpDB/Dt -evalue 1e-10 -num_threads 20 -outfmt 6 -out Dt_vs_Dt.blast"
安装MCScanX
报错
error: ‘getopt’ was not declared in this scope
dissect_multiple_alignment.cc中添加头文件
#include <getopt.h>
‘chdir’ was not declared in this scope
添加头文件 #include <unistd.h>
$ unzip MCscanX.zip
$ cd MCScanX
$ make
使用说明
http://chibba.pgml.uga.edu/mcscan2/documentation/manual.pdf
制作对应的gff文件
awk '$1~/Ghir_A/&&$3~/gene/{print $0}' Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gff3|awk -F ";" '{print $1}'|awk '{print $1"\t"substr($9,4)"\t"$4"\t"$5}'
鉴定共线性区域
-g
gap penalty
~/software/MCScanX/MCScanX -g -3 A2_A2/A2_vs_A2
下游的画点图程序
#运行dot_plotter图只能在软件的downstream那个文夹内进行
java dot_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

对duplicate gene进行分类
在输出文件中*.gene_type
概括了gff中所有基因的duplicate类型
singleton
没有基因与它重复 0proximal
in nearby chromosomal region but not adjacent 2tandem
(consecutive repeat) 串联重复 3whole genome /segmental
(i.e. collinear genes in collinear blocks) 全基因组重复 4dispersed
(other modes than segmental, tandem and proximal) 1
./duplicate_gene_classifier dir/xyz
鉴定homologous gene
1.运行All-vs-All-blastp
使用基因的第一个转录本的氨基酸序列进行blastp
在不同基因组之间进行blastp时,先将两个基因组对应的fasta文件合并为单个fasta文件,使用合并后的fasta文件分别做query和database
bsub -J "D5_vs_Dt" -q "smp" -n 20 -R span[hosts=1] -e %J.err -o %J.out "blastp -query ${DBDir}/D5_vs_Dt_pep.fasta -db ${DBDir}/D5_vs_Dt -evalue 1e-10 -num_threads 20 -outfmt 6 -out D5_vs_Dt.blast"
2.运行MCScanx
-g -3
空位罚分-3-s 5
共线性区域最小基因对数目
~/software/MCScanX/MCScanX -g -3 -s 5 A2_A2/A2_vs_A2
3.筛选同源基因对
当有一基因与多个基因同源的时候,优先选择block更大的组合;所以对共线性文件进行格式上的处理
awk '$0~/Alignment/{split($0,a," ");b=substr(a[6],3)}$1~/^[^#]/{print $0"\t"b}' A2_vs_At.collinearity
输出格式如下
第一列为block ID
第二列为block里对应的gene id
最后一列为block内含有的基因对的数目
0- 0: evm.model.Ga01G1089 evm.model.Ga02G0007 2e-134 15
0- 1: evm.model.Ga01G1104 evm.model.Ga02G0010 0 15
0- 2: evm.model.Ga01G1106 evm.model.Ga02G0012 0 15
最终一行命令解决
先提取每个block的大小
整理格式
按照基因名、block大小进行逆序排序
uniq忽略前3列进行去重
两列基因都需要进行去重
awk '$0~/Alignment/{split($0,a," ");b=substr(a[6],3)}$1~/^[^#]/{print $0"\t"b}' A2_vs_D5.collinearity |awk -F "\t" '$2~/evm/&&$3~/Gor/{split($1,a,"-");print a[1]"\t"$2"\t"$3"\t"$5}$2~/Gor/&&$3~/evm/{split($1,a,"-");print a[1]"\t"$3"\t"$2"\t"$5}'|sort -t $'\t' -k2,2 -k4,4nr|awk -F "\t" '{print $1,$3,$4,$2}' OFS="\t"|uniq -f3|awk -F "\t" '{print $1,$3,$4,$2}' OFS="\t"|sort -k4,4 -k2,2nr -t $'\t'|uniq -f3|awk '{print "Aligment_"$1"\tBlockCount_"$2"\t"$3"\t"$4}'|wc -l
4.筛选四个基因组都同源的基因
A2_vs_At_vs_D5
根据A2与At的同源关系找到对应的与A2同源的D5
A2_vs_D5_vs_Dt
根据D5与Dt的同源关系找到对应的与D5同源的A2
把两个文件A2_vs_At_vs_D5、A2_vs_D5_vs_Dt 中共有的基因对找到
根据共有的A2 vs D5 基因对找到对应于四倍体的At、Dt
将找到的At、Dt基因对与A2_vs_D进行比较,过滤一遍
使用python脚本提取四组同源基因
python ~/scripte/Alternative/module/homolog/identifyFourHomolog.py -A2D5 A2_vs_D5_collinearity.txt -AtDt At_vs_Dt_collinerity.txt -A2At A2_vs_At_collinearity.txt -D5Dt D5_vs_Dt_collinearity.txt -o A2_D5_At_Dt_collinearity.txt
Last updated
Was this helpful?