66 lines
3.7 KiB
Plaintext
Executable File
66 lines
3.7 KiB
Plaintext
Executable File
# 顺式作用元件预测，启动子预测
|
||
## 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 scripts/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 scripts/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
|
||
|
||
|
||
|