關于基因家族鑒定及表達分析步驟届惋，包括數(shù)據(jù)的獲取郑藏，軟件的使用击碗，步驟流程。需要批量的重復步驟可使用腳本來完成。

1.1 此頁面下載 hmm文件

image.png

1.2 blastp 比對鑒定putative gene

1.2.1 在擬南芥數(shù)據(jù)庫下載蛋白序列,并獲得 ubox的蛋白id埂奈，使用seqkit 或 seqtk 提取序列芹敌，不能根據(jù)序列id 提取所有序列，因為序列文件中的locus 全為大寫，而 ubox 蛋白id 個別為大寫捆等。

blast+ 可直接下載断部，解壓后為可執(zhí)行文件达址，建庫疆虚，并搜索。得到的序列較多锄贷。對石榴所有蛋白序列建庫柔昼，其他物種此家族蛋白序列做query

1.3 smart 搜索 結(jié)果序列結(jié)構(gòu)域確定目標基因 及其他特征的分析

直接使用hmmsearch的結(jié)果 用 tbtools 的smart插件做檢測蹦魔，96條序列有93條都含有此結(jié)構(gòu)域

再用pfam網(wǎng)頁測試下都含有 此結(jié)構(gòu)域。使用 SMART的結(jié)果獲得序列

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}'

根據(jù) gff文件獲得所有蛋白id與基因id的對應關系 并獲得基因id 的位置信息

grep -v "#" GCF_007655135.1_ASM765513v2_genomic.gff |awk -F "[\t=;,:]" 'BEGIN{OFS="\t"}$3=="CDS"{print$15,$17}'

有的基因有多個轉(zhuǎn)錄本讳推，怎么選取一個轉(zhuǎn)錄本

使用excel 根據(jù)基因列刪除重復項，擴展所選區(qū)域即可

將gene id 與蛋白id 染色體信息整合到一張表

根據(jù)gene在染色體位置命名

BUSCA網(wǎng)站預測亞細胞定位洒忧，ProtParam預測理化性質(zhì)（使用之前的腳本）蛉抓，驗證一條序列正確.
并將基因的一些信息整合到一起

根據(jù)gene在染色體上的信息畫圖 一般在線網(wǎng)站就可以，下載svg結(jié)果文件可使用AI編輯

1.4 MEME 分析保守基序 expasy 批量預測蛋白理化性質(zhì)

MEME suit 網(wǎng)站預測motif 一般設置motif 數(shù)量參數(shù) 及 site distribution參數(shù)（一般選anr）芝雪。在結(jié)果頁面可下載xml文件通過tbtools可提取,motif位置信息

批量預測蛋白理化性質(zhì)

1.5 比對及建樹

（1） 使用MAFFT比較準確减余，MAFFT有在線和本地版，linux 下用conda 安裝conda install -c bioconda mafft惩系，使用在線版比對

(2)有教程說用Gblocks 提取比對后的保守區(qū)域，但此方法文件引用較少如筛，就用全序列比對吧

（3）ModeFinder 計算最優(yōu)模型堡牡，其文獻中說 ModelFinder is implemented in IQ-TREE version 1.5.4 (http://www.iqtree.org) ，此外也有Modetest 杨刨。 iqtree為可執(zhí)行文件.

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搜索,注意比對結(jié)果的格式妖胀，要是后續(xù)使用  RA × ML ‐NG應使用 RA × ML ‐NG 支持的比對格式進行最優(yōu)模型查找
本來使用clustal的格式進行最佳模型查找芥颈，但ra ml ng不支持clustal格式進行建樹，又重新使用fasta格式進行模型查找赚抡，兩者所查找到的最優(yōu)模型一致 LG+F+R5爬坑。
本次實驗的 命令 nohup iqtree -s 56renameproteins.fasta -m MF -mtree -T 25 &

iqtree 后綴文件里會有最佳模型

最佳模型 LG+F+R5

（4）使用 RA × ML ‐NG 建樹，下載可執(zhí)行文件即可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}

為官方建議一般參數(shù)選擇

首先 檢查輸入文件準確性

/root/app/ramlng/raxml-ng --parse --msa example.phy --model TIM2+F+I+G4
--parse 會適合 large alignment ,且會壓縮輸入文件并給出建議的線程及內(nèi)存,結(jié)果如下盾计，如建樹指定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 數(shù)量
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

查看主機線程數(shù)量

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 &

建樹比模型查找需要較少的時間。

建樹結(jié)果文件

本次實驗的命令為

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個小時完成建樹，可能只因為參數(shù)不對總共跑了20個樹樣本

結(jié)果日志文件里說是

Starting ML tree search with 20 distinct starting trees
但是應該是1000條的哭尝，這個起始樹是啥意思哥攘，查了下代碼應該沒錯
可能需要加 --bootstrap參數(shù)

又重新 加 --bootstrap參數(shù) ，命令為

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 沒有最優(yōu)樹的結(jié)果

使用--all 參數(shù)測試 弟翘，感覺就是在試錯. 應加all參數(shù)運行

nohup /root/raxml-ng --all --msa /root/treegte/56renameproteins.fasta --model LG+F+R5 --threads 8 --bs-trees 1000 &

使用-all參數(shù)得到樹文件含有 bootstrap support values, 有很多個結(jié)果文件，可在官網(wǎng) 查看結(jié)果文件的含義莽红。其中$PREFIX.raxml.support Best-scoring ML tree with bootstrap support values 文件是最終的樹文件妥畏。

（5）itol 美化
計劃，根據(jù)不結(jié)構(gòu)域?qū)π蛄蟹纸M安吁。
給不同的基因設置不同的背景醉蚁，參照 模板，主要是設置分隔符鬼店，基因列网棍，type(標簽顏色，分枝顏色等),顏色列（顏色可從提供色板的網(wǎng)頁獲得）

1.6 protein structure

之前通過SMART獲得的domain很多妇智，先寫各腳本提取需要的domain,并重新命名序列滥玷。還真是夠猶豫的，不知道要選擇那些結(jié)果巍棱，就選數(shù)量多的惑畴，以前文獻有報道的吧。
將motif 與 domain 放在一張圖分析

1.7 順式作用原件

需要 的是基因起始上游 DNA序列航徙。以gff文件中的cds start 位置作為起始位置（測試單獨提取某個基因cds起始上游如贷，忽略 堿基 0 1 的定位問題。與此方法結(jié)果相同到踏，）杠袱。具體方法參見根據(jù)gff文件提取特定基因組序列（測試單獨提取某個基因cds起始上游，忽略 堿基 0 1 的定位問題窝稿。與此方法結(jié)果相同楣富，），其中 Feature ID只是一個序列命名的依據(jù)讹躯，搞了半天菩彬，序列長度不是2000,原來是忘記勾選某個選項

記得勾選

單個序列的html文件結(jié)果

tab文件結(jié)果

GSDS網(wǎng)站好像掛了险领，先把數(shù)據(jù)處理好吧侨舆。數(shù)據(jù)包括，cis-element的起始終止位置及不同基因上游 長度绢陌。

gsds 需要2個輸入文件挨下，一個是畫基因的骨架文件。一個是fea
用GSDS畫cis element 也看不出差別脐湾，畫成熱圖展示數(shù)量

1.8 家族內(nèi)成員的共線性

synteny 與 collonearity區(qū)別

MCScanx與 Circos 什么關系臭笆？,好像 circos是可視化 mcscanx的結(jié)果
MCScanX做共線性分析用法 中說，blast用蛋白序列或cds序列都可以

MCScanX reads in two data files: xyz.blast and xyz.gff秤掌，分別通過blastp 和從原gff文件獲得
（1） MCScanx安裝
為make后的可執(zhí)行文件
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，會生成以下文件

MCScanx 編譯后生成的文件

 /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 比對
用單物種的所有蛋白序列机杜，比對剛剛構(gòu)建的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椒拗，但起始終止位置是基因的起始終止位置。
因為比對的時候是用的所有蛋白序列获黔，里面有某個基因的轉(zhuǎn)錄本蚀苛，是不是應該取某個基因的單個轉(zhuǎn)錄本，然后再做blastp玷氏。有12831條蛋白序列是與其他蛋白序列同屬于某些基因的轉(zhuǎn)錄本堵未。這樣會減輕服務器工作量。在某個基因id只保留一個轉(zhuǎn)錄本時盏触，也要保證保留了家族中選擇的轉(zhuǎn)錄本渗蟹。
只前從gff文件中提取geneid與proteinid的對應關系時块饺，因為第二列為genome與為reseq（為葉綠體）的geneid與proteinid 位置不同。在保留單個基因的單個轉(zhuǎn)錄本時提取的序列的數(shù)量與id數(shù)量不一樣雌芽，就不用葉綠體蛋白序列做后續(xù)分析樂授艰。
并重新建庫并比對。下班世落。

后續(xù)處理得到gff文件淮腾，可根據(jù)gene的行提取，且只提取染色體上的對應gene后續(xù)將基因id替換為 單個基因id 保留的單個蛋白id.
(5) mcscan的結(jié)果處理
collinearity后綴應該含有的是共線行基因?qū)π畔⑻爰眩捎谥敖◣毂葘r的蛋白序列含有非染色體上的蛋白谷朝，現(xiàn)在只保留染色體上的蛋白序列之前的共線信息
獲得所有在染色體上的蛋白id，將染色體上非家族內(nèi)的具有共線特征的基因?qū)μ崛〕鰜恚ɑ驅(qū)Σ皇?個都是某個家族）及家族內(nèi)的共線性對提取出來
之后用Tbtools 共線性教程提供的方法畫圖武花，需要準備a,染色體長度信息圆凰，b,共線性gene對的染色體位置信息，c,一些feature的位置信息体箕，比如標注家族內(nèi)基因名稱专钉。
之后根據(jù)前提取的家族內(nèi)外的基因?qū)Γ╩cscanx有這個功能detect_collinearity_within_gene_families.pl ，與手動查找的一樣干旁。驶沼。≌海快多了）回怜，
獲得家族內(nèi)外的link region信息，并在家族 內(nèi)的link region信息后添加顏色信息. 最后使用Tbtools畫圖.

MCScanX_h與MCScanX相似不過其輸入文件不是BLASTp的結(jié)果换薄，是其他第三方軟件的結(jié)果玉雾。

按照教程即可，注意染色體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

后續(xù)寫家族內(nèi)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） 可實現(xiàn)批量對具有復制情況的基因?qū)τ嬎鉱aks, 程序中調(diào)用比對程序（需另下載）及kaks_calcultor(需另下載)复旬，不需要一對一對地比對計算了。具體操作參照paraat批量計算kaks
parat會將每對id 的kaks單獨生成一個文件冲泥。如下

parat結(jié)果

1.9 使用轉(zhuǎn)錄組數(shù)據(jù)對gene進行定量

（1）獲得基因家族的transcriptome
網(wǎng)上說是cdna序列，但數(shù)據(jù)只有 rna_from_genomic.fna 這種數(shù)據(jù)闽巩，就只能钧舌，先獲得蛋白id與rna id的對應關系担汤。
根據(jù)gff文件及rna序列文件獲得重新命名后的rna（cdna?）序列

（2）rna seq數(shù)據(jù)下載
sratoolkits 下載后為可執(zhí)行文件，不過有一個圖形選擇界面的Configuration步驟
使用sratoolkit里的prefetch 指定包含sra id的文件下載數(shù)據(jù)洼冻，fasterq-dump 將srq解壓崭歧。可shell腳本批量解壓,具體參見批量fasterqdump
服務器壞了碘赖，打算先買個按量計費的服務器驾荣。

nohup /root/sratoolkit.2.10.8-centos_linux64/bin/prefetch -p --option-file ./20sra_id.txt &
下載 需要的sra文件，

使用以下shell 腳本批量轉(zhuǎn)換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 &

最后結(jié)果都輸出到-O 指定的文件夾中普泡， 19個sra數(shù)據(jù)解壓后播掷，有418G，500gb存儲夠sra文件及轉(zhuǎn)換后文件可能不夠
（3）使用kallisto 定量

需要用到 kallisto 軟件
可直接使用conda 安裝,其他方式參見官網(wǎng)

conda install kallisto

kallisto使用,
主要用到kallisto的index 與quant命令撼班。使用腳本批量quant
Building an index: 使用基因家族的cDNA 創(chuàng)建索引歧匈，
Quantification ：根據(jù)索引，和 rna-seq數(shù)據(jù) 定量砰嘁。
rnaseq數(shù)據(jù)可為gz格式 或 plain-text

A 對輸入文件建index 索引

-i 參數(shù)指定 索引名稱
kallisto index -i 56gene_rna.idx 56gene_rna.fna

B quant 定量件炉，
需要寫個for loop, 讀取所有的sra,并并添加路徑得到每個fastq的變量.
由于quant的結(jié)果文件為以下,結(jié)果文件為統(tǒng)一的固定文件名文件，所以不能-o 指定結(jié)果輸出到一個文件夾矮湘，不然會覆蓋

quant結(jié)果

更改命令重新測試的時候斟冕，上次的目錄還是有結(jié)果文件生成，原來是直接用kill pid 只能殺死 子進程缅阳。需要先kil ppid 再kill pid磕蛇，第三次測試，重新學習標準輸出十办，標準輸入知識 數(shù)據(jù)重定向

最終使用的批量做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 然后 運行秀撇，可得到結(jié)果

nohup bash quant.sh &

按量計費，大概花費40元向族。不過要是白天連續(xù)操作應該會少點呵燕，主要是下載數(shù)據(jù)流量費用。

（4）quant結(jié)果處理
quant的結(jié)果為以下件相，即每個測序數(shù)據(jù)的quant結(jié)果保存在以sra id 命名的文件夾中再扭，tsv格式文件里有需要的信息

quant結(jié)果

由于每個sra數(shù)據(jù)的結(jié)果為單獨的一個文件夾，文件里也沒有具體的sra id 標注夜矗，需要將不同文件集合在同一文件霍衫，并不同sra結(jié)果文件里標注 sra id。就是獲得所有結(jié)果文件夾路徑侯养，進入路徑，打印tsv文件中的某些列澄干，并打印這個sra 結(jié)果文件的路徑逛揩，追加到某個文件即可柠傍。知識點，在awk 中打印 shell變量辩稽。

#進入kallisto的包含不同測序數(shù)據(jù)的結(jié)果文件夾中
for i in `find  -name "SRR*"`
do
cd $i
#echo $i #打印單個結(jié)果文件路徑
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

由于數(shù)據(jù)是寬數(shù)據(jù)？需要轉(zhuǎn)換成長數(shù)據(jù)逞泄？
接著轉(zhuǎn)換allsra.tpm文件得到患整，特定熱圖軟件的輸入格式
以tbtools輸入格式為例，得到固定輸入格式
需要將kallisto quant格式

kallisto quan結(jié)果格式

轉(zhuǎn)換為sra id 在第一行的格式喷众，即

tbtools 熱圖工具輸入格式

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)#字符串不用引號表示

就可得到輸入格式文件到千。
將不同類型轉(zhuǎn)錄組數(shù)據(jù)分開畫圖昌渤，因為處理不一樣，還是要linux 與win的行尾符不一樣憔四。還是先把結(jié)果得到吧膀息，怎么寫再看情況寫吧。
行標準化了赵，列標準化參見tbtools中標準化問題潜支，行標準化后，可以看出某基因柿汛，在不同樣本中的表達情況冗酿，這時某個樣本中的不同基因表達情況是無法比較的。列標準化苛茂，可比較已烤，某個樣本中不同目標基因的表達情況