精华内容
下载资源
问答
  • 基因家族分析

    2021-06-17 20:41:27
    #下载拟南芥基因组信息 #wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz #wget ftp://ftp.ensemblgenomes.org/pub/...
    
    #下载拟南芥基因组信息
    #wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
    #wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
    #wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
    #wget ftp://ftp.ensemblgenomes.org/pub/plants/release-41/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.41.gff3.gz
    #
    #解压压缩文件
    #gunzip *gz
    
    #观察数据,ID保持一致,也就是gff中第9列,ID标签和parent标签与蛋白序列和cds序列里面的ID是否一致;
    #处理GFF 文件里面ID中一些不必要的信息,gene:  transcript: 删除;与蛋白质中的ID保持一致:Arabidopsis_thaliana.TAIR10.pep.all.fa 
    
    #sed 's#gene:##' Arabidopsis_thaliana.TAIR10.41.gff3|sed  's#transcript:##' |sed 's#CDS:##' >Arabidopsis_thaliana.TAIR10.41.gff31
    #mv Arabidopsis_thaliana.TAIR10.41.gff31 Arabidopsis_thaliana.TAIR10.41.gff3
    
    #获取基因与mRNA的对应关系,后面会用到;
    perl script/mRNAid_to_geneid.pl Arabidopsis_thaliana.TAIR10.41.gff3 mRNA2geneID.txt
    perl script/geneid_to_mRNAid.pl Arabidopsis_thaliana.TAIR10.41.gff3 geneid2mRNAid.txt
    
    #第一次搜索结构域
    hmmsearch --domtblout WRKY_hmm_out.txt --cut_tc WRKY.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa
    
    #提取结构域序列,最后的evalue值根据实际情况可调,注意脚本提取的是第一个domain,如要提取其他domain,请修改脚本27行$a[9]==1为第一个,$a[9]==2为第二个,依次类推
    perl script/domain_xulie.pl WRKY_hmm_out.txt Arabidopsis_thaliana.TAIR10.pep.all.fa WRKY_domain.fa 1.2e-28
    
    ###########以下部分为建立物种特异模型再次搜索,可根据自己基因家族情况选做这部分内容#############################
    #clusterW多序列比对快捷方法
    echo "1\nWRKY_domain.fa\n2\n1\nWRKY_domain.aln\nWRKY_domain.dnd\nX\n\n\nX\n" |clustalw
    
    #利用比对结果建立物种特异hmm模型
    hmmbuild WRKY_domain_new.hmm WRKY_domain.aln
    
    #新建物种特异hmm模型,再次搜索
    
    hmmsearch --domtblout WRKY_domain_new_out.txt --cut_tc WRKY_domain_new.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa
    
    ############################################################################################################
    
    #筛选 hmm搜索结果,可以用excel手动筛选,筛选标准,
    #筛选EValue  <0.001
    #如果只想用hmmer搜索一次,可将下面的文件:WRKY_domain_new_out.txt 替换成第一次hmmer搜索生成的文件:WRKY_hmm_out.txt
    grep -v "^#" WRKY_domain_new_out.txt|awk '$7<0.001 {print}' >WRKY_domain_new_out_selected.txt
    
    
    #去除重复的hmmer搜索的转录本ID,多个转录本ID保留一个作为基因的代表,此步建议对脚本输出的文件手动筛选,挑选ID:
    perl script/select_redundant_mRNA.pl mRNA2geneID.txt WRKY_domain_new_out_selected.txt WRKY_remove_redundant_IDlist.txt
    
    #请手动挑选完mRNA的ID放在第一列,也就是挑选一个转录本ID代表这个基因,存成新的文件WRKY_removed_redundant_IDlist.txt:
    #利用脚本得到对应基因的蛋白序列:
    perl script/get_fa_by_id.pl WRKY_removed_redundant_IDlist.txt Arabidopsis_thaliana.TAIR10.pep.all.fa WRKY_pep_need_to_confirm.fa
    
    
    #将上面WRKY_pep_need_to_confirm.fa文件中的蛋白序列,再手动验证一下,把不需要的ID删除,最终确认的ID存成新文件:WRKY_removed_redundant_and_confirmed_IDlist.txt
    #手动确认结构域,CDD,SMART,PFAM
    #确定分子量大小:http://web.expasy.org/protparam/
    #perl script/stat_protein_fa.pl WRKY_pep_need_to_confirm.fa WRKY_pep_need_to_confirm.MW.txt
    #三大数据库网站,筛选之后去除一些不确定的基因ID,最终得到可靠的基因家族基因列表,存储在文件:WRKY_removed_redundant_and_confirmed_IDlist.txt ; 
    
    #脚本提取hmm结果文件,重新筛选一下hmm的结果:
    
    #get_data_by_id.pl 会读取第一个文件的第一列ID,然后到第二个文件中去筛选对应ID(第二个文件第一列作为ID)的数据输出到第三个文件中
    perl script/get_data_by_id.pl WRKY_removed_redundant_and_confirmed_IDlist.txt WRKY_domain_new_out_selected.txt WRKY_domain_new_out_removed_redundant.txt
    
    #截取得到序列上的保守结构域序列,注意脚本提取的是第一个domain,如要提取其他domain,请修改脚本27行$a[9]==1为第一个,$a[9]==2为第二个,依次类推
    
    perl script/domain_xulie.pl WRKY_domain_new_out_removed_redundant.txt Arabidopsis_thaliana.TAIR10.pep.all.fa WRKY_domain_confirmed.fa 0.1
    
    #得到对应基因的蛋白序列全长,脚本会读取第一个文件的第一列的ID信息,根据ID提取相应的序列:
    
    perl script/get_fa_by_id.pl WRKY_removed_redundant_and_confirmed_IDlist.txt Arabidopsis_thaliana.TAIR10.pep.all.fa WRKY_pep_confirmed.fa
    
    #得到对应基因的cds序列,脚本会读取第一个文件的第一列的ID信息,根据ID提取相应的序列:
    
    perl script/get_fa_by_id.pl WRKY_removed_redundant_and_confirmed_IDlist.txt Arabidopsis_thaliana.TAIR10.cds.all.fa WRKY_cds_confirmed.fa
    
    
    
    
    
    
    ########################进化树分析##########################################
    #cd $workdir 回到工作路径
    mkdir gene_tree_analysis
    cd gene_tree_analysis
    cp ../WRKY_domain_confirmed.fa .
    cp ../WRKY_pep_confirmed.fa .
    cp ../WRKY_cds_confirmed.fa .
    cp ../WRKY_domain_new_out_removed_redundant.txt .
    
    
    #########################利用meme软件做motif分析################################33
    
    #cd $workdir
    mkdir meme_motif_analysis
    cd meme_motif_analysis
    #搜索结构域:
    #-nmotifs 10  搜索motif的总个数
    #-minw 6   motif的最短长度
    #-maxw 50   motif的最大长度
    
    /biosoft/meme/meme-v4.12.0/bin/meme ../WRKY_pep_confirmed.fa -protein -oc ./ -nostatus -time 18000 -maxsize 6000000 -mod anr -nmotifs 10 -minw 6 -maxw 100
    
    #蛋白质保守结构域搜索
    								输出结果	因马可夫模型	蛋白序列
    hmmsearch --cut_tc --domtblout HSP_xian.txt HSP20.hmm HSP_pep_confirmed.fa
    
    ##################################基因结构分析structure####################
    
    #回到工作路径 my_gene_family
    
    mkdir gene_structure_analysis
    cd gene_structure_analysis
    cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .
    
    
    
    #获得基因的在染色体上的外显子,CDS,UTR位置信息,用于绘制基因结构图
    #注意脚本读取的是第一个(-in1)文件第一列信息,里面是转录本ID;然后把转录本对应的位置cds utr等信息筛选出来
    perl ../script/get_gene_exon_from_gff.pl -in1 WRKY_removed_redundant_and_confirmed_IDlist.txt -in2 ../Arabidopsis_thaliana.TAIR10.41.gff3 -out gene_exon_info.gff
    
    
    ################################基因定位到染色体###############################################
    # 回到工作路径  my_gene_family
    mkdir map_to_chr
    cd map_to_chr
    cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .    #WRKY基因家族ID列表文件
    
    
    #获得基因的在染色体上的位置信息,用于绘制染色体位置图,注意提取的是基因位置还是mRNA位置,以下代码是提取的 mRNA位置
    perl ../script/get_gene_weizhi.pl -in1 WRKY_removed_redundant_and_confirmed_IDlist.txt -in2 ../Arabidopsis_thaliana.TAIR10.41.gff3 -out mrna_location.txt
    
    #获得基因组染色体长度:
    samtools faidx ../Arabidopsis_thaliana.TAIR10.dna.toplevel.fa
    
    cp ../Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai .
    
    #绘图参考:http://www.omicsclass.com/article/397
    
    
    
    #######################共线性分析,基因加倍与串联重复分析MCScanX##################################
    #回到工作路径  my_gene_family
    
    mkdir MCScanX
    cd MCScanX
    
    mkdir mcscan
    
    #获取基因对应的蛋白序列信息,注意:基因的第一个转录本为代表序列;
    #从geneid2mRNAid.txt文件中每个基因挑一个转录本ID存储到allmRNAID.txt(一列)
    
    #获取基因组基因对应转录本的位置信息
    perl ../script/get_mRNA_position.pl allmRNAID.txt ../Arabidopsis_thaliana.TAIR10.41.gff3 AT.gff
    
    #获取对应蛋白序列
    perl ../script/get_fa_by_id.pl allmRNAID.txt ../Arabidopsis_thaliana.TAIR10.pep.all.fa pep.fa
    
    #blast建库
    makeblastdb -in pep.fa  -dbtype prot -title pep.fa
    
    #blastall  比对,所有基因对所有基因
    blastall -i pep.fa -d pep.fa -e 1e-10  -p blastp  -b 5 -v 5 -m 8 -o mcscan/AT.blast
    cp AT.gff mcscan/AT.gff
    
    
    #运行MCSCANX  查找共线性
    /biosoft/MCScanX/MCScanX/MCScanX mcscan/AT
    
    #下载示例文件,自己分析时需要用自己研究物种的染色体文件,和前面鉴定基因家族基因列表文件
    #wget http://chibba.pgml.uga.edu/mcscan2/examples/family.ctl
    #wget http://chibba.pgml.uga.edu/mcscan2/examples/MADS_box_family.txt
    
    sed -i 's#at##g' family.ctl
    
    #基因家族共线性绘图,注意MCSCAN绘图要回到:/biosoft/MCScanX/MCScanX/downstream_analyses 才能绘图
    cd /biosoft/MCScanX/MCScanX/downstream_analyses 
    sudo java family_circle_plotter -g /home/manager/gene2019/my_gene_family/MCScanX/mcscan/AT.gff -s /home/manager/gene2019/my_gene_family/MCScanX/mcscan/AT.collinearity -c /home/manager/gene2019/my_gene_family/MCScanX/family.ctl -f /home/manager/gene2019/my_gene_family/MCScanX/MADS_box_family.txt -o /home/manager/gene2019/my_gene_family/MCScanX/WRKY.circle.PNG
    
    
    #找到我们分析的基因家族的共线性分析关系(大片段复制现象):同样要回到mcscan的安装路径,才能分析
    
    cd /biosoft/MCScanX/MCScanX/downstream_analyses 
    perl /biosoft/MCScanX/MCScanX/downstream_analyses/detect_collinearity_within_gene_families.pl -i /home/manager/gene2019/my_gene_family/MCScanX/MADS_box_family.txt -d /home/manager/gene2019/my_gene_family/MCScanX/mcscan/AT.collinearity -o /home/manager/gene2019/my_gene_family/MCScanX/WRKY.collinear.pairs
    
    #从MCSCAN分析的结果文件:AT.tandem 提取串联重复基因可以看,:http://www.omicsclass.com/article/399
    cd /home/manager/gene2019/my_gene_family/MCScanX
    perl ../script/get_tandem_gene.pl -id gene_list.txt -tandem mcscan/AT.tandem -od ./ -name WRKY
    
    ########基因加倍分析绘图,circos###############
    #cd $workdir 回到工作路径
    mkdir circos
    cd circos
    
    cp ../MCScanX/mcscan/AT.collinearity .   
    cp ../MCScanX/WRKY.collinear.pairs .
    
    #准备circos绘图数据文件,脚本从gff里面获得位置信息并整理数据
    perl ../script/colline_v3.pl -gff ../MCScanX/AT.gff -list WRKY.collinear.pairs -colline AT.collinearity -od data -name WRKY
    
    #绘图,主要准备config3.txt配置文件,以及染色体长度文件等等。
    cd data
    
    /biosoft/circos/circos-0.69-6/bin/circos -conf config3.txt -outputdir ./ -outputfile WRKY
    
    
    
    ###############################复制基因查找 blast方法 及KAKS分析#################################
    
    #回到工作路径  my_gene_family
    
    mkdir gene_duplication_kaks_blast
    cd gene_duplication_kaks_blast
    cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .
    cp ../WRKY_cds_confirmed.fa .
    #blast建库,DNA序列,all vs all 比对,结果说明见:http://www.omicsclass.com/article/505
    makeblastdb -in WRKY_cds_confirmed.fa -dbtype nucl -title WRKY_cds_confirmed.fa 
    blastall -i WRKY_cds_confirmed.fa -d WRKY_cds_confirmed.fa -p blastn -e 1e-20  -m 8 -o WRKY_cds_confirmed_blast.out
    
    #获取基因cds序列的长度:
    samtools faidx WRKY_cds_confirmed.fa
    
    perl ../script/KAKS_SHAIXUAN.pl -in1 WRKY_cds_confirmed.fa.fai -in2 WRKY_cds_confirmed_blast.out -out duplication_gene.out
    
    ###kaks分析###
    
    echo "AT1G66600.1\nAT1G66560.1" >dupid.txt
    perl ../script/get_fa_by_id.pl dupid.txt WRKY_cds_confirmed.fa dup_gene_paired1.fa
    echo "1\ndup_gene_paired1.fa\n2\n1\ndup_gene_paired1.aln\ndup_gene_paired1.dnd\nX\n\n\nX\n" |clustalw
    
    #格式转换axt  如果遇到报错not equal,可参考:http://www.omicsclass.com/article/700
    /biosoft/KaKs_Calculator2.0/src/AXTConvertor dup_gene_paired1.aln dup_gene_paired1.axt
    /biosoft/KaKs_Calculator2.0/bin/Linux/KaKs_Calculator  -i dup_gene_paired1.axt -o dup_gene_paired1.kaks.result
    
    #注意筛选基因的位置信息:在同一条染色体,且距离小于100K
    #分离时间计算:http://www.omicsclass.com/question/896
    
    ##############################不同物种之间的共线性分析#############################################
    # 回到工作路径  my_gene_family
    mkdir mcscan_between_genome
    cd mcscan_between_genome
    
    #获取基因对应的cds序列信息,注意:基因的第一个转录本为代表序列;
    #从geneid2mRNAid.txt文件中每个基因挑一个转录本ID存储到allCDSID.txt(一列)
    
    #得到基因组上所有基因的位置信息,bed文件;以及cds序列;
    perl ../script/get_mRNA_bed.pl -in1 ../Arabidopsis_thaliana.TAIR10.41.gff3 -in2  allCDSID.txt -out ATH.bed
    
    #获取基因对应的cds序列
    perl ../script/get_fa_by_id.pl allCDSID.txt ../Arabidopsis_thaliana.TAIR10.cds.all.fa ATH.cds
    
    
    #同样的道理准备,准备白菜的基因组,bed文件和,cds文件;
    
    #处理ID:
    #sed 's#gene:##' Brassica_rapa.Brapa_1.0.41.chr.gff3|sed  's#transcript:##' |sed 's#CDS:##' >Brassica_rapa.Brapa_1.0.41.chr.gff31
    #mv Brassica_rapa.Brapa_1.0.41.chr.gff31 Brassica_rapa.Brapa_1.0.41.chr.gff3
    
    perl ../script/mRNAid_to_geneid.pl Brassica_rapa.Brapa_1.0.41.chr.gff3 BrmRNA2geneID.txt
    perl ../script/geneid_to_mRNAid.pl Brassica_rapa.Brapa_1.0.41.chr.gff3 Brgeneid2mRNAid.txt
    
    
    #获取白菜基因对应的cds序列信息,注意:基因的第一个转录本为代表序列;
    #从Brgeneid2mRNAid.txt文件中每个基因挑一个转录本ID存储到BrallCDSID.txt(一列)
    
    #得到白菜基因组上所有基因的位置信息,bed文件;以及cds序列;
    perl ../script/get_mRNA_bed.pl -in1 Brassica_rapa.Brapa_1.0.41.chr.gff3 -in2  BrallCDSID.txt -out rapa.bed
    
    #获取基因对应的cds序列
    perl ../script/get_fa_by_id.pl BrallCDSID.txt Brassica_rapa.Brapa_1.0.cds.all.fa rapa.cds
    
    
    #注意,不知道为什么,这个python版本的mcscan如果ID后面有.1 运行会报错,所以我们把拟南芥和白菜的ID统一都去除一下.1
    
    sed -i 's#\.1##' rapa.cds
    sed -i 's#\.1##' ATH.cds
    sed -i 's#\.1##' rapa.bed
    sed -i 's#\.1##' ATH.bed
    
    /biosoft/miniconda/miniconda2/bin/python -m jcvi.compara.catalog ortholog ATH rapa --cscore=0.7
    #对共线性区域进行过滤
    /biosoft/miniconda/miniconda2/bin/python -m jcvi.compara.synteny screen --minsize=0 --minspan=30 --simple ATH.rapa.anchors   ATH.rapa.anchors.new
    #绘制共线性图片:准备两个配置文件为输入文件:
    
    /biosoft/miniconda/miniconda2/bin/python -m jcvi.graphics.karyotype  --format=pdf  --figsize=15x5 mcscan_seqid mcscan_layout
    
    
    ########################结合转录组分析基因家族成员表达量########################################
    
    #回到工作路径  my_gene_family
    
    mkdir rna_seq_heatmap
    cd rna_seq_heatmap
    
    cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .
    sed -i 's/\.1//' WRKY_removed_redundant_and_confirmed_IDlist.txt
    
    
    perl ../script/get_data_by_id.pl WRKY_removed_redundant_and_confirmed_IDlist.txt GSE121407_fpkm.txt gene_family_fpkm.txt
    
    
    ###########################################以下blast为可选分析内容########################################################################
    #blastp比对寻找基因家族成员,WRKY部分
    #NCBI上搜索WRKY蛋白序列:搜索条件:WRKY[title] NOT putative[title] AND plants[filter]
    #参考文献:https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-4955-8             
    
    #回到工作路径  my_gene_family
    mkdir blast
    cd blast
    
    #blast比对首先建库  #蛋白质序列
    makeblastdb -in WRKY_NCBI_pep.fasta -dbtype prot -title WRKY_NCBI_pep.fasta   
    #
    #blastp比对
    blastall -i ../Arabidopsis_thaliana.TAIR10.pep.all.fa -d WRKY_NCBI_pep.fasta -p blastp -e 1e-10 -b 1 -v 1 -m 8 -o ncbi_WRKY_blast.out 
    
    #获得比对上的候选基因,存储在wrky.fa文件中
    perl ../script/get_fa_by_id.pl ncbi_WRKY_blast.out ../Arabidopsis_thaliana.TAIR10.pep.all.fa wrky.fa
    
    #然后将 wrky.fa 提交到 NCBI CDD  SMART 上确认,把不存在结构域的基因ID可以删除
    
    #######################基因上游顺势作用原件分析#######################################
    
    #回到工作路径 my_gene_family
    
    mkdir gene_promoter
    cd gene_promoter
    cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .
    
    #得到基因在染色体上的位置,此脚本会把基因组所有的序列读入内存,如果基因组较大,可能因为内存不足使脚本运行不成功,可以分染色体分开分析:
    perl ../script/get_gene_weizhi.pl -in1 WRKY_removed_redundant_and_confirmed_IDlist.txt -in2 ../Arabidopsis_thaliana.TAIR10.41.gff3 -out mrna_location.txt
    #根据位置信息提取,promoter序列 1500
    perl ../script/get_promoter.pl ../Arabidopsis_thaliana.TAIR10.dna.toplevel.fa mrna_location.txt promoter.fa
    
    #生成 GSDS配置文件
    cat WRKY_removed_redundant_and_confirmed_IDlist.txt|awk 'BEGIN{OFS="\t"}{print $1,"0","1500","CDS","."}' >gene.bed
    #生成feature文件
    cat PlantCARE_9210__plantCARE/plantCARE_output_PlantCARE_9210.tab|grep "Arabidopsis"|awk -F"\t"  'BEGIN{OFS="\t"} {print $1,$4,$4+length($3),$2}'>feature.bed
    
    
    展开全文
  • 用tbtools基因家族分析《一》

    千次阅读 2020-11-05 21:29:12
    基因家族分析流程蛋白序列的blast比对候选蛋白序列的提取候选基因的保守结构域确定 此版本为无代码版本,依赖tbtools来实现,纯粹记录自己的学习经验。 蛋白序列的blast比对 选择拟南芥基因家族的蛋白序列文件 选择...


    此版本为无代码版本,依赖tbtools来实现,纯粹记录自己的学习经验。

    蛋白序列的blast比对

    1. 选择拟南芥基因家族的蛋白序列文件
    2. 选择你想要查找的基因家族的基因组蛋白文件
    3. 选择你要输出的文件名
    4. 选择输出格式,默认为xml,选择"Table“格式,结果更直观
    5. 点击”start“,开始运行。
      在这里插入图片描述

    候选蛋白序列的提取

    1. 对于我们获得的结果文件rice.txt,,我们选择用excel打开,对其进行排序和去重处理,得到候选的基因ID;
      在这里插入图片描述
    2. 选择TBtools的”Fasta Extract or Filter"功能来提取上述候选基因的蛋白序列,
      在这里插入图片描述
    3. 第一个红框为水稻基因组的蛋白序列,第二个为输出的文件名,第三个就是我们的上一步得到的候选基因ID;
      在这里插入图片描述

    候选基因的保守结构域确定

    1. 利用网站NCBI CD search来预测候选基因的保守结构域;
      在这里插入图片描述
      在这里插入图片描述
    2. 根据结果来筛选候选基因,比如我们挑选的是DUF629 superfamily,其中包含DUF629 superfamily的基因则为最终的候选基因。
      在这里插入图片描述
    展开全文
  • 基因家族分析《三》

    2020-11-09 15:59:23
    基因的染色体定位分析 利用tbtools的 Gene Location Visualize From GTF/GFF工具,对染色体进行位置的绘制,这种是棒状的染色体图。 若想要圈图的展示形式,可以选择Circle Gene View。

    基因的染色体定位分析

    1. 利用tbtools的 Gene Location Visualize From GTF/GFF工具,对染色体进行位置的绘制,这种是棒状的染色体图。在这里插入图片描述
      在这里插入图片描述
    2. 若想要圈图的展示形式,可以选择Circle Gene View在这里插入图片描述
      在这里插入图片描述
    展开全文
  • 谈论到直系同源基因分析的时候,大部分教程都是介绍OrthoMCL,这是2003年发表的一个工具,目前的引用次数已经达到了3000多,但这个软件似乎在2013年之后就不在更新,而且安装时还需要用到MySQL(GitHub上有人尝试从...

    谈论到直系同源基因分析的时候,大部分教程都是介绍OrthoMCL,这是2003年发表的一个工具,目前的引用次数已经达到了3000多,但这个软件似乎在2013年之后就不在更新,而且安装时还需要用到MySQL(GitHub上有人尝试从MySQL转到sqlite)。

    而OrthoFinder则是2015年出现的软件,目前已有400多引用。该软件持续更新,安装更加友好,因此我决定使用它来做直系同源基因的相关分析。

    OrthoFinder能做什么?

    OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy提到,它的优点就是比其他的直系同源基因组的推断软件准确,并且速度还快。

    此外他还能分析所提供物种的系统发育树,将基因树中的基因重复事件映射到物种树的分支上,还提供了一些比较基因组学中的统计结果。

    OrthoFinder的分析过程

    OrthoFinder的分析过程分为如下几步:

    1. BLAST all-vs-all搜索。使用BLASTP以evalue=10e-3进行搜索,寻找潜在的同源基因。(除了BLAST, 还可以选择DIAMOND和MMSeq2)
    2. 基于基因长度和系统发育距离对BLAST bit得分进行标准化。
    3. 使用RBNHs确定同源组序列性相似度的阈值
    4. 构建直系同源组图(orthogroup graph),用作MCL的输入
    5. 使用MCL对基因进行聚类,划分直系同源组

    分析流程1 OrthoFinder2在OrthoFinder的基础上增加了物种系统发育树的构建,流程如下
    1. 为每个直系同源组构建基因系统发育树
    2. 使用STAG算法从无根基因树上构建无根物种树
    3. 使用STRIDE算法构建有根物种树
    4. 有根物种树进一步辅助构建有根基因树

    基于Duplication-Loss-Coalescent 模型,有根基因树可以用来推断物种形成和基因复制事件,最后记录在统计信息中。

    分析流程2

    软件使用

    在解压缩的OrthoFinder文件目录下(安装见最后)有一个 ExampleData, 里面就是用于测试的数据集。

    orthofinder -f ExampleData -S mmseqs
    # -f 指定文件夹
    # -S 指定序列搜索程序,有blast, mmseqs, blast_gz, diamond可用

    OrthoFinder的基本使用就是如此简单,而且最终效果也基本符合需求。

    如果你想根据多序列联配(MSA)结果按照极大似然法构建系统发育树,那么你需要加上-M msa。这样结果会更加准确,但是代价就是运行时间会更久,这是因为OrthoFinder要做10,000 - 20,000个基因树的推断。

    OrthoFinder默认用mafft进行多序列联配,用fasttree进行进化树推断。多序列联配软件还支持muscle, 进化树推断软件还支持iqtree, raxml-ng, raxml。例如参数可以设置为-M msa -A mafft -T raxml.

    并行化参数: -t参数指定序列搜索时的线程数,-a指的是序列搜索后分析的CPU数。

    软件细节

    OrthoFinder提供了config.json可以调整不同软件的参数,如下是BLASTP。

    BLASTP参数

    OrthoFinder默认使用DendroBLAST发育树,也就是根据序列相似度推断进化关系。这是作者推荐的方法,在损失部分准确性的前提下提高了运算效率。当然你可以用-M msa从多序列比对的基础上进行基因树构建。如果你先用了默认的DendroBLAST,想测试下传统的MSA方法,那么也不需要重头运行,因为有一个-b参数可以在复用之前的比对结果。

    在物种发育树的推断上,OrthoFinder使用STAG算法,利用所有进行构建系统发育树,而非单拷贝基因。此外当使用MSA方法进行系统发育树推断时,OrthoFinder为了保证有足够多的基因(大于100)用于分析,除了使用单拷贝基因外,还会挑选大部分是单拷贝基因的直系同源组。这些直系同源组的基因前后相连,用空缺字符表示缺失的基因,如果某一列存在多余50%的空缺字符,那么该列被剔除。最后基于用户指定的建树软件进行系统发育树构建。结果在"WorkingDirectory/SpeciesTree_unrooted.txt"

    使用STRIDE算法从无根树中推断出有根树, 结果就是"SpeciesTree_rooted.txt".

    结果文件

    运行结束后,会在ExampleData里多出一个文件夹,Results_Feb14, 其中Feb14是我运行的日期

    直系同源组相关结果文件,将不同的直系同源基因进行分组

    • Orthogroups.csv:用制表符分隔的文件,每一行是直系同源基因组对应的基因。
    • Orthogroups.txt: 类似于Orthogroups.csv,只不过是OrhtoMCL的输出格式
    • Orthogroups_UnassignedGenes.csv: 格式同Orthogroups.csv,只不过是物种特异性的基因
    • Orthogroups.GeneCount.csv:格式同Orthogroups.csv, 只不过不再是基因名信息,而是以基因数。

    直系同源相关文件,分析每个直系同源基因组里的直系同源基因之间关系,结果会在Orthologues_Feb14文件夹下,其中Feb14是日期

    • Gene_Trees: 每个直系同源基因基因组里的基因树
    • Recon_Gene_Trees:使用OrthoFinder duplication-loss coalescent 模型进行发育树推断
    • Potential_Rooted_Species_Trees: 可能的有根物种树
    • SpeciesTree_rooted.txt: 从所有包含STAG支持的直系同源组推断的STAG物种树
    • SpeciesTree_rooted_node_labels.txt: 同上,只不过多了一个标签信息,用于解释基因重复数据。

    比较基因组学的相关结果文件:

    • Orthogroups_SpeciesOverlaps.csv: 不同物种间的同源基因的交集
    • SingleCopyOrthogroups.txt: 单基因拷贝组的编号
    • Statistics_Overall.csv:总体统计信息
    • Statistics_PerSpecies.csv:分物种统计信息

    STAG是一种从所有基因推测物种树的算法,不同于使用单拷贝的直系同源基因进行进化树构建。

    一些重要概念:

    • Species-specific orthogroup: 一个仅来源于一个物种的直系同源组
    • Single-copy orthogroup: 在直系同源组中,每个物种里面只有一个基因。我们会用单拷贝直系同源组里的基因推断物种树以及其他数据分析。
    • Unassigned gene: 无法和其他基因进行聚类的基因。
    • G50和O50,指的是当你直系同源组按照基因数从大到小进行排列,然后累加,当加入某个组后,累计基因数大于50%的总基因数,那么所需要的直系同源组的数目就是O50,该组的基因树就是G50.

    Orthogroups, Orthologs 和 Paralogs 这三个概念推荐看图理解。

    概念辨析

    如何安装?

    最快的方法

    OrthoFinder可以通过conda安装,建议为它新建一个虚拟环境

    conda create -n orthofinder orthofinder=2.2.7

    如果你愿意折腾

    你先得安装它的三个依赖工具: MCL, FastME, DIAMOND/MMseqs2/BLAST+

    MCL有两种安装方式,最简单的就是用sudo apt-get install mcl, 但是对于大部分人可能没有root权限,因此这里用源代码编译。http://micans.org/mcl/

    wget https://www.micans.org/mcl/src/mcl-latest.tar.gz
    tar xf mcl-latest.tar.gz
    cd mcl-14.137
     ./configure --prefix=~/opt/biosoft/mcl-14.137
    make -j 20 && make install

    之后是MMseqs2, 一个蛋白搜索和聚类工具集,相关文章发表在NBT, NC上。GitHub地址为https://github.com/soedinglab/MMseqs2

    wget https://github.com/soedinglab/MMseqs2/releases/download/3-be8f6/MMseqs2-Linux-AVX2.tar.gz
    tar xzf MMseqs2-Linux-AVX2.tar.gz
    mv mmseqs2 ~/opt/biosoft/

    最后安装FastME, 这是一个基于距离的系统发育树推断软件。在http://www.atgc-montpellier.fr/fastme/binaries.php下载,上传到服务器

    下载
    tar xf fastme-2.1.5.tar.gz
    cd fastme-2.1.5
    ./configure --prefix=/opt/biosoft/fastme-2.1.5
    make && make install

    BLAST+可装可不装,推荐阅读这或许是我写的最全的BLAST教程

    以上软件安装之后,都需要将其添加到环境变量中,才能被OrthoFinder调用。

    之后在https://github.com/davidemms/OrthoFinder/releases 寻找最近的稳定版本下载到本地,例如OrthoFinder v2.2.7

    tar xzf OrthoFinder-2.2.7.tar.gz
    OrthoFinder-2.2.7/orthofinder -h
    展开全文
  • 基因家族分析《二》

    2020-11-06 11:35:04
    基因家族motif分析 打开meme网页,按下图填写参数即可; [可选] 如果你想使用linux版本的meme,参数设置如下即可。 等待结果,选择MEME XML output,右键链接另存为保存即可; 3.meme结果可视化:利用tbtools工具的...
  • 以搜索ZEB的人类的基因家族成员为例 一. 搜索已知数据库 1.在NCBI Gene上搜索ZEB1  点击summary 点击domain 点击Gene就可以看到有相同结构的Gene了   2.在UniProtKB上下载protein family list ...
  • 转载:http://www.realbio.cn/news/124.html https://blog.csdn.net/seallama/article/details/43820763 ... 1. 数据库的配置 OrthoMCL的分析需要先行建立mysql账户并建立相应的数据库。关于...
  • 使用方法参考: 基于全基因组的基因家族分析(3):SlNRAMP家族基因CDS和Genomic DNA序列获取 Gene序列要获得的ID号不仅仅是seqid,而且还需要在染色体上的位置信息——起始和终止位置,以及染色体编号。这里就需要...
  • 基因家族鉴定分析实战操作手册

    千次阅读 2019-09-14 15:18:35
    基因家族鉴定分析操作手册: 基因家族 基因家族鉴定 基因家族鉴定分析总结 1.下载基因组信息文件,gff,cds,pep,fasta文件: Ensembl下载地址:http://plants.ensembl.org/index.html 下载方法可参考:...
  • 棉花XTH基因家族全基因组鉴定及进化分析,安亚茹,杨君,本研究从亚洲棉、雷蒙德氏棉和陆地棉基因组中分别鉴定出37、39和72个XTH家族基因,编码木葡聚糖内转糖苷酶/水解酶(xyloglucan endotransgly
  • miRNA基因簇和基因家族的全局分析揭示了动态和协调的表达。
  • 大豆NAC基因家族生物信息学分析,王洋,柏锡,NAC转录因子是植物特有的、具有多种生物功能的一类重要转录因子,广泛参与植物生长发育以及生物与非生物逆境应答等。本文利用生物
  • 玉米PLDs基因家族生物信息学分析,韩雨,徐晶宇,磷脂酶D(PLD)为植物中重要的一种水解酶,PLD在逆境环境下不仅参与细胞膜的老化过程和种子萌发的过程,还作为一个信号物质响应脱落酸
  • CAFE:基因家族进化的计算分析
  • 番茄ANT基因家族的全基因组鉴定、结构特征及表达分析,杨蔚然,先志强,ANT(AINTEGUMENTA)基因家族属于植物AP2/ERF这个庞大的转录因子家族中的AP2亚家族,广泛参与植物器官发育调控和生长素信号转导相关途径并发
  • 大熊猫bHLH 基因家族的识别与分析,徐俊丽,李玉芝,碱性/螺旋-环-螺旋( basic /Helix-Loop-Helix,bHLH) 蛋白质家族是数量最多、最复杂的转录因子家族之一,几乎存在于所有真核生物体内,其功�
  • 短管赤眼蜂Argonaute蛋白基因家族鉴定及表达分析,张海燕,杨志强,动物和植物Argonaute蛋白在RNA诱导的基因沉默中发挥重要作用。本研究采用生物信息学方法在全基因组水平鉴定短管赤眼蜂Argonaute蛋白基因
  • 水稻泛素结合酶基因家族的功能分析,刘鑫,张恒,植物泛素/蛋白酶体系统在植物的生长发育、形态建成和抗病反应等过程中起着重要的作用。近年的研究又表明,某些病原菌能够模拟寄�
  • 目的研究毛果杨组蛋白去乙酰化酶(HDAC)基因家族序列及其对脱落酸(ABA)的应答反应,为HDAC 家族基因的功能研究奠定基础。方法采用生物信息学方法对毛果杨HDAC家族的16种蛋白进行系统分析;采用实时荧光定量 PCR方法,...
  • 枳BOR基因家族的克隆与缺硼胁迫下的表达分析,颜廷帅,罗庆,为了揭示柑橘砧木枳[Poncirus trifoliata (L.) Raf.]BOR基因家族各成员与硼运输的关系,本实验克隆获得枳BOR基因家族成员进行生物信息学分析
  • 森林草莓CENH3基因家族的生物信息学分析,马跃,张天龙,CENH3(centromere-specific histone H3 variant)是着丝粒特异性组蛋白H3变异体。在拟南芥中CENH3基因的改造植株与正常野生型杂交能够获得单倍体�
  • 牡丹CTR1基因家族基因片段及其cDNA 3′末端的克隆与序列分析,高娟,董丽,以牡丹品种‘洛阳红’( Paeonia suffruticosa‘Luo Yang Hong’) 花瓣为材料, 用CTAB法提取总RNA,根据已报道的CTR1(Constitutive Triple ...
  • 黄萎病菌诱导下海岛棉PR蛋白基因家族的表达分析,张桂荣,张艳,PR蛋白基因家族是一类广泛参与各种生物和非生物胁迫下植物防御应答反应的病程相关蛋白。为研究PR蛋白在棉花抗黄萎病中的功能,本�
  • 全基因组序列分析显示,杨树中的扩展蛋白基因家族包含36个成员,分属于EXPA、EXPB、EXLB和EXLA 4个亚家族。系统生物信息学分析表明,杨树扩展蛋白基因具有保守的序列特征和稳定的蛋白结构,包括9种不同的内含子起源...
  • 玉米大斑病菌漆酶基因家族的克隆及差异表达分析,刘宁,马双新,【目的】通过克隆玉米大斑病菌(Setosphaeria turcica)漆酶基因,分析并检测其在生长发育及侵染过程重的表达水平,探讨玉米大斑病菌漆�
  • Hmmer软件安装 Hmmersearch使用 Ensembl和NCBI基因组下载 samtools安装 sudo apt install samtools 此shell脚本可用hmmsearch批量处理文件 需要更名,如拟南芥pep文件更名为ath.pep。 需要将pep文件和hmm文件放在...
  • print "\n" : chomp' ${i}.cdsseq | tail -n +2 > ${i}.cdstxt echo " ${i}.cds file processing are completed and generate ${i}.cdsseq." done done 用samtools的时候如果基因号在fai文件中没有找到对应的会报错...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 283
精华内容 113
关键字:

基因家族分析