基因家族docker镜像怎么使用

这篇文章主要讲解了"基因家族docker镜像怎么使用"，文中的讲解内容简单清晰，易于学习与理解，下面请大家跟着小编的思路慢慢深入，一起来研究和学习"基因家族docker镜像怎么使用"吧！

基因家族docker 镜像使用方法

##################################################进入docker虚拟机，注意自己安装的docker版本##################################################下载基因家族分析镜像#docker pull /gene-family:v1.0#docker desktop 进入虚拟机命令#docker run -m 3G --cpus 2 --rm -v D:/gene-family:/work -it /gene-family:v1.0#docker toolbox 进入虚拟机命令#docker run -m 3G --cpus 2 --rm -v /d/gene-family:/work -it /gene-family:v1.0##############################################################################数据准备#########################################################################下载拟南芥基因组信息#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.txtperl script/geneid_to_mRNAid.pl Arabidopsis_thaliana.TAIR10.41.gff3 geneid2mRNAid.txt###################################################################################333#第一次搜索结构域####################################################################################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 -e '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#############################################################################################################筛选EValue  <0.001 hmm搜索结果，可以用excel手动筛选，筛选标准，#如果只想用hmmer搜索一次，可将下面的文件：WRKY_domain_new_out.txt 替换成第一次hmmer搜索生成的文件：WRKY_hmm_out.txtgrep -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的结果,提取结构域序列，蛋白全长，cds全长，用于后续分析##################################################################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_analysiscd gene_tree_analysiscp ../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#回到工作路径 my_gene_familymkdir meme_motif_analysiscd meme_motif_analysis#搜索结构域：#-nmotifs 10  搜索motif的总个数#-minw 6   motif的最短长度#-maxw 50   motif的最大长度meme ../WRKY_pep_confirmed.fa -protein -oc ./ -nostatus -time 18000 -maxsize 6000000 -mod anr -nmotifs 10 -minw 6 -maxw 100##################################基因结构分析structure#####################回到工作路径 my_gene_familymkdir gene_structure_analysiscd gene_structure_analysiscp ../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_familymkdir map_to_chrcd map_to_chrcp ../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.facp ../Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai .#绘图参考：http://www..com/article/397#######################基因上游顺势作用原件分析########################################回到工作路径 my_gene_familymkdir gene_promotercd gene_promotercp ../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序列 1500perl ../script/get_promoter.pl ../Arabidopsis_thaliana.TAIR10.dna.toplevel.fa mrna_location.txt promoter.fa#######################共线性分析，基因加倍与串联重复分析MCScanX###################################回到工作路径  my_gene_familymkdir MCScanXcd MCScanXmkdir 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.blastcp AT.gff mcscan/AT.gff#运行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.txtsed -i 's#at##g' family.ctl#基因家族共线性绘图，注意MCSCAN绘图要回到：/biosoft/MCScanX/MCScanX/downstream_analyses 才能绘图cd /biosoft/MCScanX/MCScanX/downstream_analyses java family_circle_plotter -g /work/my_gene_family/MCScanX/mcscan/AT.gff -s /work/my_gene_family/MCScanX/mcscan/AT.collinearity -c /work/my_gene_family/MCScanX/family.ctl -f /work/my_gene_family/MCScanX/MADS_box_family.txt -o /work/my_gene_family/MCScanX/WRKY.circle.PNG#找到我们分析的基因家族的共线性分析关系（大片段复制现象）：cd /work/my_gene_family/MCScanXperl /biosoft/MCScanX/MCScanX/downstream_analyses/detect_collinearity_within_gene_families.pl -i MADS_box_family.txt -d mcscan/AT.collinearity -o WRKY.collinear.pairs#从MCSCAN分析的结果文件：AT.tandem 提取串联重复基因可以看，：http://www..com/article/399perl ../script/get_tandem_gene.pl -id gene_list.txt -tandem mcscan/AT.tandem -od ./ -name WRKY########基因加倍分析绘图，circos################cd $workdir 回到工作路径mkdir circoscd circoscp ../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 datacircos -conf config3.txt -outputdir ./ -outputfile WRKY##############################复制基因kaks分析################################################################回到工作路径  my_gene_familymkdir gene_duplication_kakscd gene_duplication_kakscp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .cp ../WRKY_cds_confirmed.fa .echo -e "AT1G66600.1\nAT1G66560.1" >dupid.txtperl ../script/get_fa_by_id.pl dupid.txt WRKY_cds_confirmed.fa dup_gene_paired1.faecho -e "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..com/article/700AXTConvertor dup_gene_paired1.aln dup_gene_paired1.axtKaKs_Calculator  -i dup_gene_paired1.axt -o dup_gene_paired1.kaks.result#分离时间计算：http://www..com/question/896##############################不同物种之间的共线性分析############################################## 回到工作路径  my_gene_familymkdir mcscan_between_genomecd 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.gff3perl ../script/mRNAid_to_geneid.pl Brassica_rapa.Brapa_1.0.41.chr.gff3 BrmRNA2geneID.txtperl ../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.cdssed -i 's#\.1##' ATH.cdssed -i 's#\.1##' rapa.bedsed -i 's#\.1##' ATH.bedpython -m jcvi.compara.catalog ortholog ATH rapa --cscore=0.7#对共线性区域进行过滤python -m jcvi.compara.synteny screen --minsize=0 --minspan=30 --simple ATH.rapa.anchors   ATH.rapa.anchors.new#绘制共线性图片：准备两个配置文件为输入文件：python -m jcvi.graphics.karyotype  --format=pdf  --figsize=15x5 mcscan_seqid mcscan_layout########################结合转录组分析基因家族成员表达量绘制热图#########################################回到工作路径  my_gene_familymkdir rna_seq_heatmapcd rna_seq_heatmapcp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .sed -i 's/\.1//' WRKY_removed_redundant_and_confirmed_IDlist.txtperl ../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_familymkdir blastcd 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可以删除

感谢各位的阅读，以上就是"基因家族docker镜像怎么使用"的内容了，经过本文的学习后，相信大家对基因家族docker镜像怎么使用这一问题有了更深刻的体会，具体使用情况还需要大家实践验证。这里是，小编将为大家推送更多相关知识点的文章，欢迎关注！