进化树：MUSCLE/IQtree软件应用

MUSCLE:

muscle -in NB-ARC.domain.fasta -out NB-ARC.domain.afa -maxiters 2

-in 输入文件

-out 输出文件名，输出文件默认为 Fasta 格式

-maxiters 最大迭代次数

在处理大数据时，作者建议设置 最大迭代次数为 2，因为超过 2 次迭代后对模型的提升很小。

IQtree

IQtree 是一款用 极大似然法 快速、准确构建 无根进化树 的软件

输入：多序列比对 的结果文件

推荐使用 iteration=2

iqtree -s NB-ARC.domain.afa -nt 50

./iqtree -s AC.fas -p More10_part.txt -m TESTNEW -nt 4 -bb 1000 -alrt 1000

iqtree -s AC.fas -m TESTNEW -nt 16 -bb 10000 -alrt 1000

PCA软件:

---

def loadDataSet(fileName, delim='\t'):

    fr = open(fileName)

    stringArr = [line.strip().split(delim) for line in fr.readlines()]

    datArr = [map(float,line) for line in stringArr]

    return mat(datArr)

---

def pca(dataMat, topNfeat=9999999): #原数据 m*n

    meanVals = mean(dataMat, axis=0)                          #均值：1*n

    meanRemoved = dataMat - meanVals               #去除均值  m*n

    covMat = cov(meanRemoved, rowvar=0)                      #协方差矩阵 n*n

    eigVals,eigVects = linalg.eig(mat(covMat))   #特征矩阵  n*n

    eigValInd = argsort(eigVals)            #将特征值排序

    eigValInd = eigValInd[:-(topNfeat+1):-1]  #仅保留p个列（将topNfeat理解为p即可） 

    redEigVects = eigVects[:,eigValInd]      # 仅保留p个最大特征值对应的特征向量，按从大到小的顺序重组特征矩阵n*p

    lowDDataMat = meanRemoved * redEigVects #将数据转换到低维空间lowDDataMat： m*p

    reconMat = (lowDDataMat * redEigVects.T) + meanVals    #从压缩空间重构原数据reconMat：  m*n

    return lowDDataMat, reconMat

————————————————

版权声明：本文为CSDN博主「fanstuck」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。

原文链接：https://blog.csdn.net/master_hunter/article/details/114301045

R语言中的library()实现PCA

见下方链接：

(3条消息) PCA主成分分析实战和可视化 | 附R代码和测试数据_生信宝典的博客-CSDN博客_pca r代码

视频总结：

左上方脚本

1、导入数据：两张表格：标准化表达矩阵和样本信息表

gene_exp <- read.table(file = 'RNASeq_downstream/data/genes.TNM.EXPR.mtrix')

sample_info <_  read.table(file = 'RNASeq-downstream/data/sample_info.txt',

header=TRUE,

row.names=1)

2、安装并运行

BiocManager::install(PCAtools)#检查安装

library(PCAtools)#加载PCAtools

p <-  pca(mat=gene_exp.metdata=sample_info,removeVar=0.1)

view(p)

3、

pca_loadings <- p[["loadings"]]

#提取某一值

pca_rotated <- p[["rotated"]]

screenplot(p)#图一，每一主成分对样本解释读的差异

biplot(p,

x='pc1',

y='pca'，

shape='stage',#

colby='strain'，#color

legendPosition='top')#图例位置，图二聚类图

#图二，画聚类图

plotloading（p）#图三影响pc1的主要成分

linux系统

grep -v "#" blat.out |awk '$3>50 && $11<1e-30 {print $2}' | sort | uniq >protein_ids.txt