#################################################
#进入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可以删除