关于基因家族鉴定及表达分析步骤，包括数据的获取，软件的使用，步骤流程。需要批量的重复步骤可使用脚本来完成。
1.1 此页面下载 hmm文件
1.2 blastp 比对鉴定putative gene
1.2.1 在拟南芥数据库下载蛋白序列,并获得 ubox的蛋白id，使用seqkit 或 seqtk 提取序列，不能根据序列id 提取所有序列，因为序列文件中的locus 全为大写，而 ubox 蛋白id 个别为大写。
blast+ 可直接下载，解压后为可执行文件，建库，并搜索。得到的序列较多。对石榴所有蛋白序列建库，其他物种此家族蛋白序列做query
1.3 smart 搜索 结果序列结构域确定目标基因 及其他特征的分析
直接使用hmmsearch的结果 用 tbtools 的smart插件做检测，96条序列有93条都含有此结构域
再用pfam网页测试下都含有 此结构域。使用 SMART的结果获得序列
1.4 染色体及基因位置等信息
从gff文件获得染色体信息
less GCF_007655135.1_ASM765513v2_genomic.gff |awk '$2=="RefSeq"'|grep 'genome=chromosome'
从gff文件获得所有gene的信息
less GCF_007655135.1_ASM765513v2_genomic.gff |awk -F "[\t:;]" 'BEGIN{OFS="\t"}$3=="gene"{print $1,$2,$3,$4,$5,$6,$7,$11}'
根据 gff文件获得所有蛋白id与基因id的对应关系 并获得基因id 的位置信息
grep -v "#" GCF_007655135.1_ASM765513v2_genomic.gff |awk -F "[\t=;,:]" 'BEGIN{OFS="\t"}$3=="CDS"{print$15,$17}'
有的基因有多个转录本，怎么选取一个转录本
使用excel 根据基因列删除重复项，扩展所选区域即可
将gene id 与蛋白id 染色体信息整合到一张表
根据gene在染色体位置命名
BUSCA网站预测亚细胞定位，ProtParam预测理化性质（使用之前的脚本），验证一条序列正确.
并将基因的一些信息整合到一起
根据gene在染色体上的信息画图 一般在线网站就可以，下载svg结果文件可使用AI编辑
1.4 MEME 分析保守基序 expasy 批量预测蛋白理化性质
MEME suit 网站预测motif 一般设置motif 数量参数 及 site distribution参数（一般选anr）。在结果页面可下载xml文件通过tbtools可提取,motif位置信息
1.5 比对及建树
（1） 使用MAFFT比较准确，MAFFT有在线和本地版，linux 下用conda 安装conda install -c bioconda mafft，使用在线版比对
(2)有教程说用Gblocks 提取比对后的保守区域，但此方法文件引用较少，就用全序列比对吧
（3）ModeFinder 计算最优模型，其文献中说 ModelFinder is implemented in IQ-TREE version 1.5.4 (http://www.iqtree.org) ，此外也有Modetest 。 iqtree为可执行文件.
If you have enough computational resource, you can perform a thorough and more accurate analysis that invokes a full tree search for each model considered via the -mtree option:
iqtree -s example.phy -m MF -mtree -T AUTO # 查找所有模型中最佳模型，-T 指定线程信息， -mtree 使用所有Model搜索,注意比对结果的格式，要是后续使用 RA × ML ‐NG应使用 RA × ML ‐NG 支持的比对格式进行最优模型查找
本来使用clustal的格式进行最佳模型查找，但ra ml ng不支持clustal格式进行建树，又重新使用fasta格式进行模型查找，两者所查找到的最优模型一致 LG+F+R5。
本次实验的 命令 nohup iqtree -s 56renameproteins.fasta -m MF -mtree -T 25 &
iqtree 后缀文件里会有最佳模型
（4）使用 RA × ML ‐NG 建树，下载可执行文件即可raxml-ng git hub](https://github.com/amkozlov/raxml-ng)
RAxML-NG supports alignments in FASTA, PHYLIP and CATG formats.，可指定格式
raxml-ng --msa prim.phy --model GTR+G --prefix T4 --threads 2 --seed 2 --tree pars{25},rand{25}
为官方建议一般参数选择
首先 检查输入文件准确性
/root/app/ramlng/raxml-ng --parse --msa example.phy --model TIM2+F+I+G4
--parse 会适合 large alignment ,且会压缩输入文件并给出建议的线程及内存,结果如下，如建树指定threads过多。会说
WARNING: You might be using too many threads (20) for your alignment with 2460 unique patterns.
NOTE: For the optimal throughput, please consider using fewer threads
所以需要 使用 /raxml-ng --parse 得到软件建议的threads 数量
RAxML-NG v. 0.9.0 released on 20.05.2019 by The Exelixis Lab.
Developed by: Alexey M. Kozlov and Alexandros Stamatakis.
Contributors: Diego Darriba, Tomas Flouri, Benoit Morel, Sarah Lutteropp, Ben Bettisworth.
Latest version: https://github.com/amkozlov/raxml-ng
Questions/problems/suggestions? Please visit: https://groups.google.com/forum/#!forum/raxml
RAxML-NG was called at 07-Aug-2020 09:46:57 as follows:
/home/Pomgroup/gdp/app/raxml_ng/raxml-ng --parse --msa 56renameproteins.fasta --model LG+F+R5
Analysis options:
run mode: Alignment parsing and compression
start tree(s):
random seed: 1596764817
tip-inner: OFF
pattern compression: ON
per-rate scalers: OFF
site repeats: ON
branch lengths: proportional (ML estimate, algorithm: NR-FAST)
SIMD kernels: AVX2
parallelization: PTHREADS (16 threads), thread pinning: OFF
[00:00:00] Reading alignment from file: 56renameproteins.fasta
[00:00:00] Loaded alignment with 56 taxa and 3929 sites
Alignment comprises 1 partitions and 2460 patterns
Partition 0: noname
Model: LG+FC+R5
Alignment sites / patterns: 3929 / 2460
Gaps: 82.41 %
Invariant sites: 42.84 %
NOTE: Binary MSA file created: 56renameproteins.fasta.raxml.rba
* Estimated memory requirements : 207 MB
* Recommended number of threads / MPI processes: 8
Please note that numbers given above are rough estimates only.
Actual memory consumption and parallel performance on your system may differ!
Alignment can be successfully read by RAxML-NG.
Execution log saved to: /home/Pomgroup/gdp/gf/56renameproteins.fasta.raxml.log
Analysis started: 07-Aug-2020 09:46:57 / finished: 07-Aug-2020 09:46:57
Elapsed time: 0.021 seconds
查看主机线程数量
grep 'processor' /proc/cpuinfo | sort -u | wc -l
建树 指定线程及 bootstrap
nohup /root/app/ramlng/raxml-ng --msa example.phy --model TIM2+F+I+G4 -threads 1 --bs-trees 1000 &
建树比模型查找需要较少的时间。
本次实验的命令为
nohup /home/Pomgroup/gdp/app/raxml_ng/raxml-ng --msa ./56renameproteins.fasta --model LG+F+R5 --threads 8 --bs-trees 1000 &
uj说可能需要1-2天。因为指定1000次重复。
Analysis started: 07-Aug-2020 09:53:15 / finished: 07-Aug-2020 10:59:47
56条序列大概1个小时完成建树，可能只因为参数不对总共跑了20个树样本
结果日志文件里说是
Starting ML tree search with 20 distinct starting trees
但是应该是1000条的，这个起始树是啥意思，查了下代码应该没错
可能需要加 --bootstrap参数
又重新 加 --bootstrap参数 ，命令为
nohup /home/Pomgroup/gdp/app/raxml_ng/raxml-ng --bootstrap --msa ./56renameproteins.fasta --model LG+F+R5 --threads 10 --bs-trees 1000 &
测试了一次，已超过24小时（27个小时可能），只跑了590多个树，56条序列也可能需要2天，但是加--bootstrap 没有最优树的结果
使用--all 参数测试 ，感觉就是在试错. 应加all参数运行
nohup /root/raxml-ng --all --msa /root/treegte/56renameproteins.fasta --model LG+F+R5 --threads 8 --bs-trees 1000 &
使用-all参数得到树文件含有 bootstrap support values, 有很多个结果文件，可在官网 查看结果文件的含义。其中$PREFIX.raxml.support Best-scoring ML tree with bootstrap support values 文件是最终的树文件。
（5）itol 美化
计划，根据不结构域对序列分组。
给不同的基因设置不同的背景，参照 模板，主要是设置分隔符，基因列，type(标签颜色，分枝颜色等),颜色列（颜色可从提供色板的网页获得）
1.6 protein structure
之前通过SMART获得的domain很多，先写各脚本提取需要的domain,并重新命名序列。还真是够犹豫的，不知道要选择那些结果，就选数量多的，以前文献有报道的吧。
将motif 与 domain 放在一张图分析
1.7 顺式作用原件
需要 的是基因起始上游 DNA序列。以gff文件中的cds start 位置作为起始位置（测试单独提取某个基因cds起始上游，忽略 碱基 0 1 的定位问题。与此方法结果相同，）。具体方法参见根据gff文件提取特定基因组序列（测试单独提取某个基因cds起始上游，忽略 碱基 0 1 的定位问题。与此方法结果相同，），其中 Feature ID只是一个序列命名的依据，搞了半天，序列长度不是2000,原来是忘记勾选某个选项
Plantcare 会发送压缩包结果至提供的邮箱，解压后有一个tab文件，及所有序列的html结果
下图为某序列的html结果的某个cis element 的信息
下图为tab文件中某序列的html结果的某个cis element 的信息
根据2个文件，tab文件中的列可能分别为序列id. cis名称， cis的保守序列，cis的起始位置，cis-element长度，正负链，搜索的数据库物种，功能。
根据 cis的名称提取所需的行。并根据功能筛选。
newplace也可预测cis element 此网站会有对cis element 的解释.
GSDS网站好像挂了，先把数据处理好吧。数据包括，cis-element的起始终止位置及不同基因上游 长度。
gsds 需要2个输入文件，一个是画基因的骨架文件。一个是fea
用GSDS画cis element 也看不出差别，画成热图展示数量
1.8 家族内成员的共线性
MCScanx与 Circos 什么关系？,好像 circos是可视化 mcscanx的结果
MCScanX做共线性分析用法 中说，blast用蛋白序列或cds序列都可以
MCScanX reads in two data files: xyz.blast and xyz.gff，分别通过blastp 和从原gff文件获得
（1） MCScanx安装
为make后的可执行文件
MCScanx安装 教程
参照 更改文件信息，给
msa.h
dissect_multiple_alignment.h
detect_collinear_tandem_arrays.h
这三个文件前面添加上#include <unistd.h>
make报错 javac: Command not found 安装yum install -y java-1.8.0-openjdk-devel.x86_64
make成功lt，会生成以下文件
（2）使用makebalstdb 对单物种所有蛋白序列建库
首先删除序列Id后的注释信息 sed命令即可
/root/app/blast/ncbi-blast-2.10.1+/bin/makeblastdb -in onepid_forallgene.faa -dbtype prot -out all -logfile allpep.log -title all
(3) blastp 比对
用单物种的所有蛋白序列，比对刚刚构建的database
此处更改e 值 及 num_alignments，记得输出格式为 6
nohup /root/app/blast/ncbi-blast-2.10.1+/bin/blastp -query onepid_forallgene.faa -db all -evalue 1e-10 -num_alignments 5 -outfmt 6 -out pg_pg.blast &
（4）gene 的gff文件
因为比对的时候是用的蛋白序列，gff中也应该是蛋白序列id，但起始终止位置是基因的起始终止位置。
因为比对的时候是用的所有蛋白序列，里面有某个基因的转录本，是不是应该取某个基因的单个转录本，然后再做blastp。有12831条蛋白序列是与其他蛋白序列同属于某些基因的转录本。这样会减轻服务器工作量。在某个基因id只保留一个转录本时，也要保证保留了家族中选择的转录本。
只前从gff文件中提取geneid与proteinid的对应关系时，因为第二列为genome与为reseq（为叶绿体）的geneid与proteinid 位置不同。在保留单个基因的单个转录本时提取的序列的数量与id数量不一样，就不用叶绿体蛋白序列做后续分析乐。
并重新建库并比对。下班。
后续处理得到gff文件，可根据gene的行提取，且只提取染色体上的对应gene后续将基因id替换为 单个基因id 保留的单个蛋白id.
(5) mcscan的结果处理
collinearity后缀应该含有的是共线行基因对信息，由于之前建库比对时的蛋白序列含有非染色体上的蛋白，现在只保留染色体上的蛋白序列之前的共线信息
获得所有在染色体上的蛋白id，将染色体上非家族内的具有共线特征的基因对提取出来（基因对不是2个都是某个家族）及家族内的共线性对提取出来
之后用Tbtools 共线性教程提供的方法画图，需要准备a,染色体长度信息，b,共线性gene对的染色体位置信息，c,一些feature的位置信息，比如标注家族内基因名称。
之后根据前提取的家族内外的基因对（mcscanx有这个功能detect_collinearity_within_gene_families.pl ，与手动查找的一样。。。快多了），
获得家族内外的link region信息，并在家族 内的link region信息后添加颜色信息. 最后使用Tbtools画图.
MCScanX_h与MCScanX相似不过其输入文件不是BLASTp的结果，是其他第三方软件的结果。
按照教程即可，注意染色体name要一致
(6) 基因再染色体上的位置
# less -SN 56gene_protein_info.txt|awk '{print $8}'|sort|uniq -c
4 chr1
5 chr2
8 chr3
8 chr4
10 chr5
10 chr6
7 chr7
4 chr8
1 chr_stop
后续写家族内collinerity 分析。 及不同物种间的情况 circos学习 也可用pyhton的一个模块画 https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)
（7）物种间共线性
# nohup /root/app/blast/ncbi-blast-2.10.1+/bin/blastp -query 64at_ubox.fa -db pg56 -evalue 1e-5 -outfmt 6 -out atpg2.blast &
(8) 计算kaks
[paraat](https://bigd.big.ac.cn/tools/paraat） 可实现批量对具有复制情况的基因对计算kaks, 程序中调用比对程序（需另下载）及kaks_calcultor(需另下载)，不需要一对一对地比对计算了。具体操作参照paraat批量计算kaks
parat会将每对id 的kaks单独生成一个文件。如下
物种间的kaks，拟南芥的蛋白序列及cds序列在拟南芥数据库网站下载。
Paraat的输入文件为 具有复制关系的基因对，所含的cds序列，及蛋白序列，其他参数看官网介绍，其中-p 指定一个文件，文件内容为指定的线程数
1.9 使用转录组数据对gene进行定量
（1）获得基因家族的transcriptome
网上说是cdna序列，但数据只有 rna_from_genomic.fna 这种数据，就只能，先获得蛋白id与rna id的对应关系。
根据gff文件及rna序列文件获得重新命名后的rna（cdna?）序列
（2）rna seq数据下载
sratoolkits 下载后为可执行文件，不过有一个图形选择界面的Configuration步骤
使用sratoolkit里的prefetch 指定包含sra id的文件下载数据，fasterq-dump 将srq解压。可shell脚本批量解压,具体参见批量fasterqdump
服务器坏了，打算先买个按量计费的服务器。
nohup /root/sratoolkit.2.10.8-centos_linux64/bin/prefetch -p --option-file ./20sra_id.txt &
下载 需要的sra文件，
使用以下shell 脚本批量转换sra文件
for i in `find -name "*.sra"`
do
echo $i
/root/sratoolkit.2.10.8-centos_linux64/bin/fasterq-dump -p -m 2800 -O /root/heatmap/allfasterq/ --split-3 $i
done
命名为fd.sh ,然后后台不挂断运行脚本
nohup bash fd.sh &
最后结果都输出到-O 指定的文件夹中， 19个sra数据解压后，有418G，500gb存储够sra文件及转换后文件可能不够
（3）使用kallisto 定量
需要用到 kallisto 软件
可直接使用conda 安装,其他方式参见官网
conda install kallisto
kallisto使用,
主要用到kallisto的index 与quant命令。使用脚本批量quant
Building an index: 使用基因家族的cDNA 创建索引，
Quantification ：根据索引，和 rna-seq数据 定量。
rnaseq数据可为gz格式 或 plain-text
A 对输入文件建index 索引
-i 参数指定 索引名称
kallisto index -i 56gene_rna.idx 56gene_rna.fna
B quant 定量，
需要写个for loop, 读取所有的sra,并并添加路径得到每个fastq的变量.
由于quant的结果文件为以下,结果文件为统一的固定文件名文件，所以不能-o 指定结果输出到一个文件夹，不然会覆盖
更改命令重新测试的时候，上次的目录还是有结果文件生成，原来是直接用kill pid 只能杀死 子进程。需要先kil ppid 再kill pid，第三次测试，重新学习标准输出，标准输入知识 数据重定向
最终使用的批量做quant的脚本为
sraid=$(cat /root/heatmap/20sra_id.txt)
echo $sraid
for sra in ${sraid}
do
f1="/root/heatmap/allfasterq/${sra}.sra_1.fastq"
f2="/root/heatmap/allfasterq/${sra}.sra_2.fastq"
outdir="/root/heatmap/quantout/${sra}/"
echo "quanting $sra"
echo ${f1} #测试路径是否正确
echo ${f2}
echo ${outdir} #输出文件路径
kallisto quant -i /root/heatmap/56gene_rna.idx -o ${outdir} -t 4 -b 50 ${f1} ${f2}
done
命名为quant.sh 然后 运行，可得到结果
nohup bash quant.sh &
按量计费，大概花费40元。不过要是白天连续操作应该会少点，主要是下载数据流量费用。
（4）quant结果处理
quant的结果为以下，即每个测序数据的quant结果保存在以sra id 命名的文件夹中，tsv格式文件里有需要的信息
由于每个sra数据的结果为单独的一个文件夹，文件里也没有具体的sra id 标注，需要将不同文件集合在同一文件，并不同sra结果文件里标注 sra id。就是获得所有结果文件夹路径，进入路径，打印tsv文件中的某些列，并打印这个sra 结果文件的路径，追加到某个文件即可。知识点，在awk 中打印 shell变量。
#进入kallisto的包含不同测序数据的结果文件夹中
for i in `find -name "SRR*"`
do
cd $i
#echo $i #打印单个结果文件路径
less abundance.tsv|awk -v a="$i" '{print a,$1,$5}' >>/root/project/gf/input/heatmapget/quantout/allsra.tpm
cd ..
done
将以上脚本命名为tpmget.sh，运行即可
bash tpmget.sh
由于数据是宽数据？需要转换成长数据？
接着转换allsra.tpm文件得到，特定热图软件的输入格式
以tbtools输入格式为例，得到固定输入格式
需要将kallisto quant格式
转换为sra id 在第一行的格式，即
需要在R 中使用以下SCRIPT,主要是dcast函数的使用
library(reshape2) #需要的包
setwd('C:/Users/Acer/Desktop/codee/rr/project/56genetpmformatchaange/')
b<-read.table("56gene_20sra.tpm",stringsAsFactors = F)
dim(b)
head(b)
tpm<-dcast(a,formula = V2~V1) #V2 为基因列，V1为sra id 列，
write.table(tpm,file="56gene_20sra_convert.tpm2",
row.names = FALSE, #不写行号
quote = FALSE)#字符串不用引号表示
就可得到输入格式文件。
将不同类型转录组数据分开画图，因为处理不一样，还是要linux 与win的行尾符不一样。还是先把结果得到吧，怎么写再看情况写吧。
行标准化，列标准化参见tbtools中标准化问题，行标准化后，可以看出某基因，在不同样本中的表达情况，这时某个样本中的不同基因表达情况是无法比较的。列标准化，可比较，某个样本中不同目标基因的表达情况
