前言
??做生分析有時候需要從基因組fasta文件中提取基因的序列，對于這個需求有不少現(xiàn)成的軟件可以來實(shí)現(xiàn)财著，今天來跟大家分享4款可以完成這個需求的軟件联四，廢話不多說了，直接來看看是哪4個軟件：seqkit撑教、seqtk朝墩、bedtools、gffread伟姐，這四個軟件本身功能都很豐富且實(shí)用收苏，提取序列只是軟件的功能之一。前面兩款功能很類似愤兵，這個從它們的名字也可以看出來鹿霸，功能有些冗余實(shí)際實(shí)用時選擇其中一款即可。bedtools應(yīng)該很多人都用過秆乳，這個軟件主要是用來操作bed文件杜跷，可以對兩個bed進(jìn)行交集、并集等矫夷，gffread單獨(dú)拿出來有些人可能不知道葛闷，如果說cufflinks很多人可能就知道了，gffread就是cufflinks工具集中的一個子軟件双藕，它還用一個常用的功能就是gtf淑趾、gff格式互轉(zhuǎn)。那下面我們就來看看四款軟件的具體使用方法忧陪。
- seqkit
首先來看一下seqkit使用方法扣泊，該軟件的功能很多(可使用seqkit -h來查看全部可用命令)，提取序列使用的是subseq子命令嘶摊，該命令可以根據(jù)bed延蟹、gtf、region位置信息從fasta文件中提取基因序列叶堆，并且可以擴(kuò)展兩側(cè)的序列阱飘，具體代碼如下：
Usage:
seqkit subseq [flags]
Flags:
--bed string by tab-delimited BED file
--chr strings select limited sequence with sequence IDs when using --gtf or --bed (multiple value supported, case ignored)
-d, --down-stream int down stream length
--feature strings select limited feature types (multiple value supported, case ignored, only works with GTF)
--gtf string by GTF (version 2.2) file
--gtf-tag string output this tag as sequence comment (default "gene_id")
-h, --help help for subseq
-f, --only-flank only return up/down stream sequence
-r, --region string by region. e.g 1:12 for first 12 bases, -12:-1 for last 12 bases, 13:-1 for cutting first 12 bases. type "seqkit subseq -h" for more examples
-u, --up-stream int up stream length
Global Flags:
--alphabet-guess-seq-length int length of sequence prefix of the first FASTA record based on which seqkit guesses the sequence type (0 for whole seq) (default 10000)
--id-ncbi FASTA head is NCBI-style, e.g. >gi|110645304|ref|NC_002516.2| Pseud...
--id-regexp string regular expression for parsing ID (default "^(\\S+)\\s?")
--infile-list string file of input files list (one file per line), if given, they are appended to files from cli arguments
-w, --line-width int line width when outputing FASTA format (0 for no wrap) (default 60)
-o, --out-file string out file ("-" for stdout, suffix .gz for gzipped out) (default "-")
--quiet be quiet and do not show extra information
-t, --seq-type string sequence type (dna|rna|protein|unlimit|auto) (for auto, it automatically detect by the first sequence) (default "auto")
-j, --threads int number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)
#提取序列
seqkit subseq --bed gencode_lncrna.bed -o test.fa GRCh38.primary_assembly.genome.fa
可以看來seqkit提取序列還是很方便的，而且可以分別指定上下游延申多少bp的側(cè)翼序列，或者只提取上下游序列也是可以的沥匈，還可以指定輸出文件的寬度蔗喂，更多細(xì)節(jié)功能這里就不多說了。
- seqtk
下面來看一下seqtk的功能高帖，這個軟件根據(jù)bed文件來從fasta中提取序列缰儿，從下面的代碼可以看出這個軟件的子命令參數(shù)很少，通過重定向的方式輸出到新的文件里散址，用起來很簡單乖阵。
Usage: seqtk subseq [options] <in.fa> <in.bed>|<name.list>
Options: -t TAB delimited output
-l INT sequence line length [0]
#提取序列
seqtk subseq GRCh38.primary_assembly.genome.fa gencode_lncrna.bed >test.fa
- bedtools
bedtools功能很多，非常出名的功能是用來對bed文件取交集等一系列操作预麸，那些功能這里就不多說了瞪浸，感興趣的可以自己去測試∈ζ椋可以看出該軟件可以根據(jù)bed默终、gff、vcf位置來從fasta中提取序列犁罩，下面來演示一下該軟件如何提取序列齐蔽。
Usage: bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf> -fo <fasta>
Options:
-fi Input FASTA file
-bed BED/GFF/VCF file of ranges to extract from -fi
-fo Output file (can be FASTA or TAB-delimited)
-name Use the name field for the FASTA header
-split given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
-tab Write output in TAB delimited format.
- Default is FASTA format.
-s Force strandedness. If the feature occupies the antisense,
strand, the sequence will be reverse complemented.
- By default, strand information is ignored.
-fullHeader Use full fasta header.
- By default, only the word before the first space or tab is used.
#提取序列
bedtools getfasta -fi GRCh38.primary_assembly.genome.fa -bed gencode_lncrna.bed -name -fo test.fa
該軟件提取序列時，可以根據(jù)基因鏈的信息來提取序列床估，這個功能有時候很實(shí)用含滴。同時這個軟件的速度也是非常的塊。
- gffread
gffread是cufflinks軟件集中的一個子軟件丐巫，該軟件另一個常用的功能是gtf谈况、gff格式互轉(zhuǎn)，感興趣的可以自己測試递胧。這個軟件的參數(shù)非常多碑韵，這里就不展示了。這款軟件可以根據(jù)gff缎脾、gtf位置信息來從fasta中提取序列祝闻，值得一提的是該軟件在提取序列時，可以選擇輸出的格式遗菠，如提取外顯子區(qū)的序列联喘、CDS序列、翻譯成蛋白序列辙纬，功能很是貼心豁遭。下面來展示一下如何提取。
#提取序列
gffread -g GRCh38.primary_assembly.genome.fa -w test.fa gencode.v34.primary_assembly.annotation.gtf
介紹了幾個可以根據(jù)bed文件來提取序列的軟件贺拣，那有些人可能還想知道如何得到bed文件呢蓖谢？可以根據(jù)gtf文件來生成捂蕴，下面給大家介紹一款用來根據(jù)gtf生成bed文件的軟件——bedops。這里提醒一下蜈抓，該軟件有很多種安裝方式启绰，但推薦大家使用conda安裝省事省心昂儒，其他方式基本都需要sudo權(quán)限沟使。下面來看一下這款軟件使用方法。
#安裝軟件
conda install -c conda-forge -c bioconda bedops
#gtf to bed
gtf2bed < gencode.v34.primary_assembly.annotation.gtf >test.bed
用conda安裝好bedops后渊跋，激活一下環(huán)境就可以使用了腊嗡，轉(zhuǎn)換格式用的是gtf2bed子軟件，值得注意的是該軟件的輸入輸出文件都用重定向的方式拾酝，還有一個需要注意的點(diǎn)提醒一下大家燕少，如果你的gtf文件中沒有transcript ID這個屬性會報錯，會出現(xiàn)這樣的錯誤提示“Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)” 蒿囤。此時客们，可通過下面的命令來解決問題：
#gtf文件沒有transcript ID屬性可以使用該命令
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' input.gtf | gtf2bed - > output.bed
最后
??今天給大家介紹了四款可以用來從fasta提取序列的軟件，外加一個可以將gtf轉(zhuǎn)bed文件的軟件材诽，有了這幾款軟件底挫，基本上提取序列就是很容易的事。這些軟件其他的功能也是很不錯脸侥，大家可以自己去探索吧建邓！emm，今天就分享到這里了睁枕。
