基因家族鉴定与分析基因家族鉴定与分析

关于基因家族鉴定及表达分析步骤，包括数据的获取，软件的使用，步骤流程。需要批量的重复步骤可使用脚本来完成。
1.1 此页面下载 hmm文件
文章图片
image.png 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 后缀文件里会有最佳模型

文章图片
最佳模型 LG+F+R5 （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/raxmlRAxML-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+R5Analysis 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 sitesAlignment comprises 1 partitions and 2460 patternsPartition 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: 8Please 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.logAnalysis started: 07-Aug-2020 09:46:57 / finished: 07-Aug-2020 09:46:57Elapsed 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网站预测cis-element.
Plantcare 会发送压缩包结果至提供的邮箱，解压后有一个tab文件，及所有序列的html结果
下图为某序列的html结果的某个cis element 的信息

文章图片
单个序列的html文件结果
下图为tab文件中某序列的html结果的某个cis element 的信息

文章图片
tab文件结果
根据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 家族内成员的共线性 synteny 与 collonearity区别
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
make报错 javac: Command not found 安装yum install -y java-1.8.0-openjdk-devel.x86_64
make成功lt，会生成以下文件

文章图片
MCScanx 编译后生成的文件
（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文件

文章图片
mcscan输入的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单独生成一个文件。如下

文章图片
parat结果
物种间的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 指定结果输出到一个文件夹，不然会覆盖

文章图片
quant结果更改命令重新测试的时候，上次的目录还是有结果文件生成，原来是直接用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格式文件里有需要的信息

文章图片
quant结果由于每个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格式

文章图片
kallisto quan结果格式
转换为sra id 在第一行的格式，即

文章图片
tbtools 热图工具输入格式
需要在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中标准化问题，行标准化后，可以看出某基因，在不同样本中的表达情况，这时某个样本中的不同基因表达情况是无法比较的。列标准化，可比较，某个样本中不同目标基因的表达情况