# 顺式作用元件预测，启动子预测
## 1.数据
- gff3
ln -s  /home/xia-lab/wutx/00test/geneFamily20200803/63.SPL_structure_motif/02.gene_structure/Ft.SPL.gff3

- genome 
ln -s  /home/xia-lab/wutx/00test/geneFamily20200803/61.SPL_identification/00.genome/Ft.genome.fasta

- geneID
ln -s  /home/xia-lab/wutx/00test/geneFamily20200803/63.SPL_structure_motif/02.gene_structure/Ft.SPL.geneID

## 2.获取plantCARE预测数据
- 数据处理1：gff格式转gtf格式，mrna都会变成transcript
gffread  -T  -o  Ft.SPL.gtf  Ft.SPL.gff3

- 数据处理2：提取基因上游2000bp序列
seqkit  subseq   --gtf Ft.SPL.gtf  --feature transcript   --up-stream 2000   --id-ncbi    --only-flank    --gtf-tag  transcript_id   Ft.genome.fasta  > gene.upstream.fasta

- 数据处理3：修订ID, >FtPinG0005035600.01.T01
awk '{if($1~/>/){print ">"$2 }else{print $0}}' gene.upstream.fasta > gene.upstream.changeID.fasta

- 上传PlantCARE预测，结果文件发送至邮箱, 24条蛋白质序列, 大致花费10分钟，并下载邮件附件至工作目录。正常步骤：
(1)打开网站https://bioinformatics.psb.ugent.be/webtools/plantcare/html/；
(2)点击Searche for CARE;
(3)在"email address to send the results back"中填写邮箱402468026@qq.com；
(4)在“reference name or ID for the sequence (optional)”填写项目名称，直接取用“Prefix”如苦荞的SPL基因家族分析取Prefix为“Ft.SPL”;
(5)在“please only submit fasta formated sequences as pure/simple text files, no word documents or other binary things”中点击浏览文件，选择前述数据处理3的结果文件“gene.upstream.changeID.fasta”上传;
(6)点击Search按钮提交；
(7)等到邮箱收到来自no-reply@psb.vib-ugent.be发送的邮件名为“no-reply”的附件(PlantCARE_566_spl_plantCARE.tar.gz,即*.plantCARE.tar.gz的附件即为预测结果)；

- 解压PlantCARE预测结果文件*.plantCARE.tar.gz，过滤成bed格式，plantCARE_output_PlantCARE_14593.tab.转换 plantCARE_output_PlantCARE_*.tab格式成bed格式, plantcare结果需要进行手动过滤，挑选绘图的element 
awk -F "\t" '{print $1"\t"$4-1"\t"$4+length($3)"\t"$2}' plantCARE_output_PlantCARE_14593.tab  > PlantCARE.bed

### 上传GSDS绘图
- 生成promoter序列的bed文件，用于GSDS软件展示pep的染色体bed， 基因的2000bp长度bed
ln -s   ../../02.gene_structure/Ft.SPL.geneID
awk '{print $1"\t0\t2000\tCDS\t."}' Ft.SPL.geneID >   gene.upstream.bed

- 绘图：基因gene.upstream.bed和MYC启动子的bed.由于plantCARE 预测结果会比较多，根据需要对PlantCARE.bed 进行下过滤再绘图。
grep "MYC" PlantCARE.bed > PlantCARE_MYC.bed

### R语言绘图
- 提取元件描述文件1：统计元件的频数
awk '{print $2}' plantCARE_output_PlantCARE_14593.tab | sort | uniq -c | sort -nr > element_counts.txt

- 提取元件描述文件2：统计元件在基因上的频数
awk '{count[$1"\t"$2]++} END {for(key in count) print key "\t" count[key]}' \
    plantCARE_output_PlantCARE_14593.tab > gene_element_counts.txt
    
- 提取元件描述文件3：提取元件+描述列
python3 tiqu_miaoshuwenjian.py --input cir_element.desc.yuanyou.csv

- R绘图1：数据整理
(1)树文件Ft.SPL.2.iqtree.treefile
sed 's/\01.T01//g' Ft.SPL.1.iqtree.treefile > Ft.SPL.2.iqtree.treefile
(2)亚家族分类文件Ft.SPL.geneID.lable.group.2.subfamily,只保留一串id，不需要有.01.T01之类的后缀
sed 's/\.01.T01//g' Ft.SPL.geneID.lable.group.subfamily > Ft.SPL.geneID.lable.group.2.subfamily
(3)功能文件表格
cir_element.desc.yuanyou.csv

- R绘图2：绘图
Rscript promoter_analysis_complete.R Ft.SPL.2.iqtree.treefile Ft.SPL.geneID.lable.group.2.subfamily ../plantCARE_output_PlantCARE_21318.tab cir_element.desc.yuanyou_clean.csv Ft_SPL_analysis



