格式

fasta是一种仅包含序列与序列名称的文件格式，主要用于存储生物的序列信息，如基因组、基因的核酸序列或氨基酸序列，文件拓展名一般为fa、fasta、fna，是最常见的生物数据格式。fasta的序列格式如下，每条序列的信息包括两行：

>NP_224538.1 hypothetical protein
ATGACGGGTTCTGTTGACCGGCCCGACCAGAATCGCGGTGAGCGATCAATGAAGTCACCAGGGTTGGATTTGGTCAGGCGCACCCTGGACGAAGCTCGTGCTGCTGCCCGCGCGCGCGGACAAGACGCCGGTCGAGGGCGGGTCGCTTCCGTTGCGTCGGGTCGGGTGGCCGGACGGCGACGAAGCTGGTCGGGTCCGGGGCCCGACATTCGTGATCCACAACCGCTGGGTAAGGCCGCTCGTGAGCTGGCAAAGAAACGCGGCTGGTCGGTGCGGGTCGCCGAGGGTATGGTGCTCGGCCAGTGGTCTGCGGTGGTCGGCCACCAGATCGCCGAACATGCACGCCCGACTGCGCTAAACGACGGGGTGTTGAGCGTGATTGCGGAGTCGACGGCGTGGGCGACGCAGTTGAGGATCATGCAGGCCCAGCTTCTGGCCAAGATCGCCGCAGCGGTTGGCAACGATGTGGTGCGATCGCTAAAGATCACCGGGCCGGCGGCACCATCGTGGCGCAAGGGGCCTCGCCATATTGCCGGTAGGGGTCCGCGCGACACCTACGG
ATAA

第一行：以“>”开头的任意文字说明，用于标记序列，是每条序列的唯一标志符，以保证后续分析中能区分文件中的序列，序列ID中可以包含注释信息。
第二行：使用既定的核苷酸或氨基酸编码符号代表的序列，可以在一行，也可以被分割成多行，直到下一个>开头的行之前都是一条序列。

获取

ncbi Genome，ebi上存储有大量的基因组或者基因序列，Uniport上存储有大量的蛋白组或蛋白序列，可以从这些地方下载。当获得下载链接地址后，使用wget命令直接下载（Windows随意）。

wget 超链接地址

处理

文件与序列处理

#格式化
seqtk  seq -l 0 gene.fa #所有序列显示在一行
seqtk  seq -l 50 gene.fa #每行显示50个碱基
## 排序 使用seqkit sort通过id/名称/序列/长度来排序
seqkit sort [flags] #Flags: -l,通过序列长度；-n,通过全名而不是id；-s,通过序列；-i,忽视大小写；-r,反向排序；-2,双流程模式读取文件来降低内存使用。
##统计
seqkit stats XX.fa #统计序列条数，碱基总数，reads读长分布等
seqtk comp XX.fa #统计fastq文件每条序列ATCG四种碱基组成以及质量值分布
seqtk  seq -l 0 gene.fa  | grep -v ">" | awk '{print length($0)}' | head #统计每条序列长度
seqtk  seq -l 0 gene.fa  |grep -v ">" | awk '{print length($0)}' | sort | uniq -c #统计长度并按照长度计算频数
##提取序列
  #  lcl|NC_000962.3_cds_NP_216564.2_2038
  #  lcl|NC_000962.3_cds_YP_177964.1_3341
seqtk subseq gene.fa name.list #根据ID列表提取序列
##截取序列
  #  lcl|NC_000962.3_cds_NP_216564.2_2038    100    2000
  #  lcl|NC_000962.3_cds_YP_177964.1_3341     45     567
seqtk subseq gene.fa name.list #根据包含｜ID 起始位置 终止位置｜的bed文件截取序列
##按长度过滤
seqtk seq -L 1000 gene.fa #过滤长度小于1000bp序列
seqkit seq -m 1000 gene.fa #过滤长度小于1000bp序列
seqkit seq -M 1000 gene.fa #过滤长度大于1000bp序列
##反向互补
seqtk seq -r test.fa #seqtk取反向互补序列
seqkit seq -r test.fa #seqkit取反向序列
seqkit seq -r -p test.fa #seqkit seq 加-r -p同时取反向互补序列
##转换大小写
seqkit seq -l gene.fa #小写
seqkit seq -u gene.fa #大写

功能相关处理

##预测基因
prodigal -a ref.pep -d ref.cds -f gff -g 11  -o ref.gff  -s ref.stat -i ref.fna >prodigal.log #原核生物
  #  -i：输入文件，fasta格式
  #  -o：输出结果文件，有多种格式可选
  #  -f：输出文件类型gbk, gff, or sco
  #  -d：基因的核酸序列
  #  -a：基因的氨基酸序列
  #  -g：密码子表，细菌为第11
  #  -p：模式，单菌还是宏基因组
  #  -s：统计信息
augustus --strand=both --genemodel=partial --singlestrand=false --protein=on --introns=on --start=on --stop=on --cds=on --codingseq=on --alternatives-from-evidence=true --gff3=on --UTR=on --outfile=XX.gff --species=human human.genome.fa #真核生物（案例为人类）
  #  --strand=both/forward /backward 表示注释基因在两条链还是其中一条
  #  --genemodel=partial/intronless/complete/atleastone /exactlyone  
  #  partial : 允许在序列边界预测不完整的基因(默认)   intronless : 只预测单外显子基因，适用于原核生物某些些真核生物中  complete : 只预测完整基因  atleastone : 预测至少一个完整的基因   exactlyone : 准确预测一个完整的基因
  #  --singlestrand=true 独立预测每条链上的基因，允许在相反的链上有重叠的基因，默认关闭
  #  --protein=on/off；--introns=on/off；--start=on/off；--stop=on/off；--cds=on/off；--codingseq=on/off
  #  输出选项。输出预测蛋白序列，内含子，起始密码子，终止密码子。或者在“初始”、“内部”、“终端”和“单外显子”之外使用“cds”。cds不包括停止密码子(除非stopCodonExcludedFromCDS=false)，而终端和单个外显子包含停止密码子
  #  --alternatives-from-evidence=true report alternative transcripts when they are suggested by hints
  #  --gff3=on/off 输出gff3格式
  #  --UTR=on/off 预测除编码序列外的未翻译区域。目前，这只适用于人类，galdieria, toxopl asma和caenorhabditis
  #  --outfile=filename 打印输出到文件
  #  --species 物种
##预测rRNA
rnammer -S bac -m tsu,lsu,ssu  -gff ref.gff -f ref.frn ref.fna  
  #  -S：物种类型，古细菌，细菌或者真菌
  #  -m：需要rRNA类型，如果真要16S，则单独选择lsu  
  #  -gff：输出gff格式结果
  #  -f：输出fasta格式序列
##预测tRNA
tRNAscan-SE  -B  -o tRNAScan.out -f tRNAScan.out.structure -m stat.list ref.fna  
  #  -B ：物种为细菌； -A ：物种为古细菌；-O ：输入序列为细胞器； -G ：包括全部类型
  #  -o：输出结果
  #  -f：tRNA二级结构
  #  -m：统计结果
##翻译氨基酸
有多种在线工具例如ExPaSy(https://web.expasy.org/translate/)，EMBOSS Transeq from EBI(https://www.ebi.ac.uk/Tools/st/emboss_transeq/)，DNA to ProteinTranslation(http://bio.lundberg.gu.se/edu/translat.html)等支持将核酸序列转换成蛋白质序列，直接上传基因的核酸序列即可。
##功能注释
emapper.py -m diamond -i XX.fa -o Prefix --output data_dir --cpu 6 #eggnog-mapper注释
  #  -i：输入文件，基因的氨基酸序列
  #  -m：选择运行模式
  #  -h：输出帮助文档
  #  --output：输出结果前缀
  #  --output_dir：输出结果目录
##查找重复序列
RepeatMasker -pa 2 -q -species tuber -html -gff -dir repeat ref.fna
  #  -pa：线程数
  #  -q：快速模式，敏感性稍低，-s为慢速模式，敏感性更高
  #  -species：物种名
  #  -html：输出html结果
  #  -gff：输出gff格式结果
  #  -dir：输出文件夹
##串联重复序列
trf ref.fna 2 7 7 80 10 50 500 -f -d -m 
  #  2 7 7 80 10 50：为运行模式选项，各种罚分标准。
  #  -m：输出屏蔽序列
  #  -f ：输出侧翼序列
  #  -d ：输出结果文件
  #  -h：输出html格式结果
##基因组拼接 fasta格式的pacbio或者nanopore测序，可以直接进行拼接。
canu -p ouput -d result genomeSize=2.8m -num_threads 2 -nanopore-raw pacbio.fasta 
##建立索引
bwa index -a is ref.fna #bwa建立索引

##blast比对
makeblastdb -in ref.fa -dbtype nucl -parse_seqids -out ref.fa #核酸序列建立索引
blastn -query gene.fa -db ref.fna -out blastn.out -outfmt 0 -evalue 1e-5 #blastn比对
————————————————————————————————————————————————————————————————————————————————————————
makeblastdb -in XX.fasta -dbtype prot -out XX.db #蛋白质序列建立索引
blastp -query test.fasta -out test.blast -db XX.db -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_alignments 1 -num_threads 4 #blastp比对

##多序列比对 用于构建系统发育树，对多个样品的fasta可以直接使用muscle，clustalw，mafft，megacc进行多序列比对。
muscle -in input1.fasta input2.fasta input3.fasta -out output.mfa #muscle
orthofinder -f fasta_dir/ -S diamond -M msa -T raxml -t 6 #orthofinder

———以上纯属个人理解与记录