使用Linux版的MEGA构建某一基因家族的基因进化树

最近碰到了个需求，让我构建某一基因家族的进化树，并根据进化关系进行相应的分类。这让我想起了之前上课的时候，一个做进化的老师给我们讲过，如果不是纯做进化方向的课题，MEGA完全够用了。由于windows的内存等有限，做几十个基因还凑合，要是上百个基因就吃不消了，于是就想到了用Linux下的MEGA来做。

https://www.megasoftware.net/

由于是二进制文件，直接解压缩，添加到环境变量就可以用了。

具体请看我这篇文章。 https://www.jianshu.com/p/5fd60c818651

上一步我得到了该基因家族的所有基因家族的蛋白序列，然后我用windows下的MEGA的 muscle 算法进行了比较，【align-build alignment-上一步的基因家族蛋白序列-muscle比对-data-export-FASTA format】

最终我得到了比对后的 multiproteins.fasta 文件。

首先进行参数的解读，相比于其他软件，我觉得这款软件比较好理解，也容易上手。

这里的 .mao 文件尤为重要，较为简单的方法是拿到windows下去设置，具体请看组学大讲堂的这篇推送。

https://www.omicsclass.com/article/568

版本信息

MEGA version 10.1.8

For 64-bit Linux

Build 10200331

参数解读

EXAMPLES

This example performs a multiple sequence alignment on codons (it assumes that you have created the file "Clustal_Codon_Alignment.mao" using the prototyper (megaproto). A fasta file with coding data is used as input and the resulting alignment is output in the MEGA format:

This example shows how to construct a neighbor-joining phylogeny for each of a list of sequence data files.

The analysis will be performed for each file listed in "listOfDataFiles.txt" and all results will be written to

the ~/Documents/outputDirectory/ directory:

megacc -a ~/Documents/NJ_Tree_Settings.mao -l ~/Documents/listOfDataFiles.txt -o ~/Documents/outputDirectory/

LIST FORMAT

When using the -l option, each file to be analyzed must be on its own line. For example:

~/Documents/myData/seqData1.fas

~/Documents/myData/seqData2.fas

~/Documents/myData/seqData3.fas

If the analyses are to use a user-provided Newick tree file, then the tree files are given on the same line as the data files, following two pipe characters. For example:

~/Documents/myData/seqData1.fas || ~/Documents/myData/treeFile1.nwk

~/Documents/myData/seqData2.fas || ~/Documents/myData/treeFile2.nwk

~/Documents/myData/seqData3.fas || ~/Documents/myData/treeFile3.nwk

我的最终使用：

下一步我打算用 ggtree 来美化，具体学习情况，我再更新。

Gblocks它可以消除DNA或蛋白质序列中排列不一致的位置和不同的区域。这些位置可能不是同源的，或者可能已经被多个取代所饱和，在进行系统发育分析之前消除比较方便。Gblocks选择块的方式与通常手工 *** 作的方式类似，但遵循可重复的一组条件。对于缺少大段连续的非保守位置、间隙位置缺乏或密度低、侧翼位置保守程度高的情况，选择的块体必须满足一定的要求，使最终的对齐更适合于系统发育分析。Gblocks输出几个文件以可视化所选的块。Gblocks等程序的使用减少了手工编辑多个比对的必要性，使自动化大数据集的系统发育分析成为可能，最后，便利了比对的复制和其他研究人员后续的系统发育分析。

linux使用：

Gblocks proteins.fasta -t=p   # 使用蛋白质序列 ，

# protenins.fasta 指需要处理的fasta format，也可以是NBRF/PIR format 

更多参数如下：

-t=p/d/c  # 序列类型p/d/c分别对应蛋白质、DNA、密码子

-b1=       #保守位置的最小序列数(50%的序列数+ 1)，应该大于序列数的一半，小于或等于序列总数的任意整数

-b2=       #最小数量的序列为侧翼位置，(序列数的85%)等于或大于保守位置序列的最小数目的任意整数

详情见 ：http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_documentation.html

trimmomatic是一款用来处理illumina测序数据的工具，可以是单条的single reads，也可以是成对的pairend reads。支持压缩格式数据。功能和其他数据处理的程序都差不多，主要包括，1、去除adapter序列以及测序中其他特殊序列；2、采用滑动窗口的方法，切除或者删除低质量碱基 1. 先新建一个文件夹，mkdir trimmomatic 2.  cd Trimmomatic   （后ls） 3. wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip 4. unzip Trimmomatic-0.38.zip 5. cd Trimmomatic-0.38  （后ls） 6.  which java  (java 在/opt/tesc/share/jdk1.8.0-20/bin/java中） 7.  /opt/tesc/share/jdk1.8.0-20/bin/java   （后ls） 8. pwd 9. /opt/tsce/share/jdk1.8.0_20/bin/java -jar /home/HYZ930402/Zmq/trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar 10.  ls 11. /opt/tsce/share/jdk1.8.0_20/bin/java -jar /home/HYZ930402/Zmq/trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar  --help 12. 进入自己fastq数据的文件夹 13./opt/tsce/share/jdk1.8.0_20/bin/java -jar /home/HYZ930402/Zmq/trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 017_R1.fastq 017_R2.fastq output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 （/opt/tsce/share/jdk1.8.0_20/bin/java为Java的路径，-jar /home/HYZ930402/Zmq/trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar为该软件所在的位置，需要明确指明质量值体系是Phred33还是Phred64，默认是Phred64，这需要特别注意，因为我们现在的测序数据基本都是Phred33的了，所以一定要指定这个参数。017_R1.fastq 017_R2.fastq要进行过滤的文件，output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz为输出文件 ILLUMINACLIP，接头序列切除参数。LLUMINACLIP:TruSeq3-PE.fa:2:30:10（省掉了路径）意思分别是：TruSeq3-PE.fa是接头序列，2是比对时接头序列时所允许的最大错配数；30指的是要求PE的两条read同时和PE的adapter序列比对，匹配度加起来超30%，那么就认为这对PE的read含有adapter，并在对应的位置需要进行切除【注】。10和前面的30不同，它指的是，我就什么也不管，反正只要这条read的某部分和adpater序列有超过10%的匹配率，那么就代表含有adapter了，需要进行去除； LEADING，规定read开头的碱基是否要被切除的质量阈值； TRAILING，规定read末尾的碱基是否要被切除的质量阈值； SLIDINGWINDOW，滑动窗口长度的参数，SLIDINGWINDOW:5:20代表窗口长度为5，窗口中的平均质量值至少为20，否则会开始切除； MINLEN，规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉。14. 若要将002_R1.fastq改为002_R1.fastq.gz，直接gzip 002_R2.fastq即可，若要解压，直接gunzip 002_R2.fastq.gz trimmomatic可以对测序数据进行过滤 java -jar trimmomatic-0.35.jar PE -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:$file_path/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 运行上面的命令可以完成以下任务 Remove adapters (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10) #去掉接头 Remove leading low quality or N bases (below quality 3) (LEADING:3) #去掉开头质量低于3或N碱基 Remove trailing low quality or N bases (below quality 3) (TRAILING:3) #去掉末尾质量低于3或N碱基 Scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 15 (SLIDINGWINDOW:4:15) #以4个碱基为滑动窗口对reads进行扫描，当平均质量值低于15时进行剪切 Drop reads below the 36 bases long (MINLEN:36) #去掉长度小于36的reads

---