基因家族docker 镜像使用方法

基因家族docker 镜像使用方法
#################################################
#进入docker虚拟机,注意自己安装的docker版本
#################################################
#下载基因家族分析镜像
#docker pull omicsclass/gene-family:v1.0
#docker desktop 进入虚拟机命令
#docker run -m 3G --cpus 2 --rm -v D:/gene-family:/work -it omicsclass/gene-family:v1.0
#docker toolbox 进入虚拟机命令
#docker run -m 3G --cpus 2 --rm -v /d/gene-family:/work -it omicsclass/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.txt
perl 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.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的结果,提取结构域序列,蛋白全长,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_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
#回到工作路径 my_gene_family
mkdir meme_motif_analysis
cd 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_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
#######################基因上游顺势作用原件分析#######################################
#回到工作路径 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
#######################共线性分析,基因加倍与串联重复分析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  查找共线性
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 
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/MCScanX
perl /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.omicsclass.com/article/399
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
circos -conf config3.txt -outputdir ./ -outputfile WRKY

##############################复制基因kaks分析###############################################################
#回到工作路径  my_gene_family
mkdir gene_duplication_kaks
cd gene_duplication_kaks
cp ../WRKY_removed_redundant_and_confirmed_IDlist.txt .
cp ../WRKY_cds_confirmed.fa .
echo -e "AT1G66600.1\nAT1G66560.1" >dupid.txt
perl ../script/get_fa_by_id.pl dupid.txt WRKY_cds_confirmed.fa dup_gene_paired1.fa
echo -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.omicsclass.com/article/700
AXTConvertor dup_gene_paired1.aln dup_gene_paired1.axt
KaKs_Calculator  -i dup_gene_paired1.axt -o dup_gene_paired1.kaks.result
#分离时间计算: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
python -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_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可以删除


0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

654 篇文章

作家榜 »

  1. omicsgene 654 文章
  2. 安生水 325 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. CORNERSTONE 72 文章
  7. rzx 67 文章
  8. xun 66 文章