[基因组工具]seqkit的使用 |
您所在的位置：网站首页 › 太铁佳苑地址查询 › [基因组工具]seqkit的使用 |
|
SeqKit的学习 --20191017
软件的介绍
SeqKit是一种跨平台的、极快的，全面的fasta/q处理工具。SeqKit为所有的主流操作系统提供了一种可执行的双元文件，包括Windows，Linux，Mac OS X，并且不依赖于任何的配置或预先配置就可以直接使用。 软件的安装 ## Install via conda conda install -c bioconda seqkit 软件的命令 ## 序列和子序列 **seq** 转换序列（序列颠倒，序列互补，提取ID） **subseq** 从区域/gtf/bed中获得序列，包括侧面的序列 **sliding** 滑动序列，支持环式基因组 **stats** 对FASTA/Q files进行简单统计 **faidx** 创造fasta索引文件并提取子序列 **watch** 检测并连线序列特点的柱状图 **sana** 清除质量不好的单线的fastq文件 ## 格式转换 **fx2tab** 将FASTA/Q 文件转变成表格形式 (1th: name/ID, 2nd: sequence, 3rd: quality) **tab2fx** 转变表格形式为fasta/q格式 **fq2fa** 转变fastq文件为fasta文件 **convert** 在Sanger, Solexa and Illumina中转换fastq的质量编码 **translate** 将DNA/RNA序列转变成蛋白序列（支持模棱两可的碱基） ## 搜索 **grep** 根据ID/名称/序列/序列motif 搜索序列，且允许错配 **locate** 定位子序列/motif,且允许错配 **fish** 使用本地比对在较大序列中寻找短序列 **amplicon** 经由引物检索扩增子（或它附近特定的区域） ## bam文件的处理和监视 **bam** 监视和连线bam文件记录特点的直方图 ## 设置参数 **head** 打印第一个Nfasta/q的记录 **range** 在一个范围内（start:end）打印fasta/q的记录 **sample** 通过数量或比例来体验序列 **rmdup** 通过id/名称/序列 来去除复制的序列 **duplicate** 复制N次的序列 **common** 通过id/名称/序列 发现多条序列中共有的序列 **split** 通过id/seq region/size/parts (mainly for FASTA) 将序列劈开成文件 **split2** 将序列通过大小或部分 劈开成文件 ## 编辑 **replace** 通过规律表达来代替名字或序列 **rename** 重新命名复制的ID **restart** 为环状基因组重新设置起始位置 **concat** 从多个文件中经由相同的ID来连接序列 **mutate** 编辑序列（点突，插入，删除） ## 排序 **shuffle** 变换序列位置 **sort** 将序列经由id/name/sequence 进行排序 软件命令详解Sequence ID 大部分的软件，包括seqkit默认将主导的非空格字母作为ID。 FASTA header ID >123456 gene name 123456 >longname longname >gi|110645304|ref|NC_002516.2| Pseudomona gi|110645304|ref|NC_002516.2| 举例说明软件如何使用 ##下载参考序列，一个fastq文件，两个fasta文件 wget http://data.biostarhandbook.com/reads/duplicated-reads.fq.gz wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.protein.faa.gz对fastq文件进行一个概括浏览 $ seqkit stat *.gz file format type num_seqs sum_len min_len avg_len max_len duplicated-reads.fq.gz FASTQ DNA 15,000 1,515,000 101 101 101 viral.1.1.genomic.fna.gz FASTA DNA 6 195,842 5,386 32,640.3 154,675 viral.1.protein.faa.gz FASTA Protein 112 55,855 38 498.7 3,122在fasta/q文件中获取每条序列的GC含量 GC content $ seqkit fx2tab --name --only-id --gc viral*.fna.gz NC_001798.2 70.38 NC_030692.1 50.10 NC_027892.1 40.57 NC_029642.1 39.88 NC_001607.1 50.07 NC_001422.1 44.76 Custom bases (A, C and A+C) $ seqkit fx2tab -H -n -i -B a -B c -B ac viral.1.1.genomic.fna.gz #name seq qual a c ac NC_001798.2 14.87 35.03 49.90 NC_030692.1 25.03 25.25 50.28 NC_027892.1 29.68 19.44 49.11 NC_029642.1 30.14 19.22 49.36 NC_001607.1 25.30 25.54 50.84 NC_001422.1 23.97 21.48 45.45 附seqkit fx2tab 的使用（用来统计序列的信息） Usage: seqkit fx2tab [flags] Flags: -a, --alphabet 打印字母表字母 -q, --avg-qual 打印read的平均质量 -B, --base-content strings 输出指定的碱基含量 -g, --gc 输出GC含量 -H, --header-line 打印标题列 -l, --length 打印序列的长度 -n, --name 只打印名字，而没有序列或者质量 -i, --only-id 只打印基因的ID从fastq/a中根据名字和ID提取序列子集 $ seqkit sample --proportion 0.001 duplicated-reads.fq.gz \ | seqkit seq --name --only-id > id.txt ##管道符前面的命令是随机取总文件0.1%的序列，管道符后面的是提取前面的符合要求的序列的ID。 ## 查看ID list内容 $ head id.txt SRR1972739.2996 SRR1972739.3044 SRR1972739.3562 ##通过ID list文件来搜索序列 $ seqkit grep --pattern-file id.txt duplicated-reads.fq.gz > duplicated-reads.subset.fq.gz seqkit sample 的使用 -n, --number int 通过数量随机提取序列（结果也许并不完全匹配） -p, --proportion float 通过比例随机提取序列 -s, --rand-seed int 随机函数 -2, --two-pass 减少内存占用 举例：随机抽取序列 seqkit sample -n 10000 -s 11 test1_1.fq -o sample.fq seqkit sample -p 0.1 -s 11 test1_1.fq -o sample.fq从fasta/q序列中找到合并碱基并找到它的位置(这个仿佛有点难度，不报错也不打印内容到屏幕)？？ $ seqkit fx2tab -n -i -a viral.1.1.genomic.fna.gz \ | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]" ## 定位简并碱基，e.g, N and K，这个仿佛也有问题，不报错也不打印内容到屏幕？？？ seqkit grep --pattern-file id2.txt viral.1.1.genomic.fna.gz \ | seqkit locate --ignore-case --only-positive-strand --pattern K+ --pattern N+移去相同序列中重复的fasta/q记录 $ seqkit rmdup --by-seq --ignore-case duplicated-reads.fq.gz > duplicated-reads.uniq.fq.gz [INFO] 7172 duplicated records removed seqkit rmdup 的使用（用来通过id/名称/序列来去除重复的序列） Usage: seqkit rmdup [flags] Flags: -n, --by-name 通过全名而不是id -s, --by-seq 通过序列 -D, --dup-num-file string 保存数量的文件并列出重复的seqs -d, --dup-seqs-file string 保存重复seqs的文件 -i, --ignore-case 忽视字母大小写在fastq/a序列中定位motif/子序列/酶切位点 $ cat enzymes.fa >EcoRI GAATTC >MmeI TCCRAC >SacI GAGCTC >XcmI CCANNNNNNNNNTGG $ seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz 输出的内容 seqID patternName pattern strand start end matched gi|526245010|ref|NC_021865.1| MmeI TCCRAC + 1816 1821 TCCGAC gi|526245010|ref|NC_021865.1| SacI GAGCTC + 19506 19511 GAGCTC gi|526245010|ref|NC_021865.1| XcmI CCANNNNNNNNNTGG + 2221 2235 CCATATTTAGTGTGG seqkit locate 的使用 Usage: seqkit locate [flags] -d, --degenerate 包含简并碱基模式和motif --gtf 输出为GTF格式 -i, --ignore-case 忽视字母大小写 -m, --max-mismatch int 通过序列匹配时允许的最大错配 -G, --non-greedy 非贪婪模式，更快但是可能错过与其他重叠的motif -P, --only-positive-strand 只搜索正链 -f, --pattern-file 模式或motif文件（fasta格式） -p, --pattern strings 搜索pattern或motif seqkit locate -i -d -p AUGGACUN test.fa 输出结果： seqID patternName pattern strand start end matched cel-mir-58a AUGGACUN AUGGACUN + 81 88 AUGGACUG ath-MIR163 AUGGACUN AUGGACUN - 122 129 AUGGACUC怎样通过长度来对大量的fasta文件进行排序 $ seqkit sort --by-length viral.1.1.genomic.fna.gz > viral.1.1.genomic.sorted.fa seqkit sort 的使用 通过id/名称/序列/长度来排序，对于fasta格式的文件，使用flag -2 来减少内存的使用，不支持fastq格式的文件 Usage:seqkit sort [flags] Flags: -l, --by-length 通过序列长度 -n, --by-name 通过全名而不是id -s, --by-seq 通过序列 -i, --ignore-case 忽视大小写 -r, --reverse 反向排序 -2, --two-pass 双流程模式读取文件来降低内存使用根据标题信息来拆分fasta序列 ## 先对fasta文件的标题进行概括浏览 $ seqkit head -n 3 viral.1.protein.faa.gz | seqkit seq --name --only-id YP_009137150.1 YP_009137151.1 YP_009137152.1 ## 将默认的ID转换成新的ID $ seqkit head -n 3 viral.1.protein.faa.gz \ | seqkit seq --name --only-id --id-regexp "\[(.+)\]" 这一步并不懂？？？ Human alphaherpesvirus 2 Human alphaherpesvirus 2 Human alphaherpesvirus 2 ## 根据header进行拆分,会生成一个文件夹 $ seqkit split --by-id --id-regexp "\[(.+)\]" viral.1.protein.faa.gz seqkit head 的使用（首先打印第N条fasta/q记录） Usage:seqkit head [flags] Flags: -n, --number int 先打印N个fasta/q的记录（默认是10） seqkit seq 的使用 对序列进行转换（颠倒，互补，提取ID等） Usage:seqkit seq [flags] Flags: -p, --complement 取互补序列 --dna2rna DNA到RNA --rna2dna RNA到DNA -G, --gap-letters string gap letters (default "- \t.") -l, --lower-case 用小写字母打印序列 -M, --max-len int 只打印短于最大长度的的序列 -n, --name 只打印name -g, --remove-gaps 移去组装序列中的gap -r, --reverse 取反向序列 -i, --only-id 只打印ID而不是全名 -q, --qual 只打印品质？？？ -s, --seq 只打印序列 --id-regexp string 对于解析ID的正则表达式，(default "^(\\S+)\\s?") seqkit split 的使用 usage：seqkit split [flags] -i, --by-id 根据序列ID切割序列 -p, --by-part int 将一份文件分割成N份 -s, --by-size int 将一个文件按照N条序列进行分割 -O, --out-dir string 输出路径 -2, --two-pass 降低内存的使用从一个text文件已知的字符串中搜索并替换fasta标题 $ seqkit grep --by-name --use-regexp --ignore-case --pattern hypothetical viral.1.protein.faa.gz \ | seqkit head -n 3 | seqkit seq --name ## 此命令也是有报错，说明没有执行好？？？ $ seqkit replace --kv-file changes.tsv --pattern "^([^ ]+ )(\w+) " \ --replacement "\${1}{kv} " --key-capt-idx 2 --keep-key viral.1.protein.faa.gz > renamed.fa seqkit grep 的使用（通过ID/名称/序列/序列motif来搜索序列，允许错配） Usage:seqkit grep [flags] Examples: 1-based index 1 2 3 4 5 6 7 8 9 10 negative index 0-9-8-7-6-5-4-3-2-1 seq A C G T N a c g t n 1:1 A 2:4 C G T -4:-2 c g t -4:-1 c g t n -1:-1 n 2:-2 C G T N a c g t 1:-1 A C G T N a c g t n 1:12 A C G T N a c g t n -12:-1 A C G T N a c g t n -n, --by-name 通过全名来匹配而不是ID -s, --by-seq 搜索seq中的子seq，只搜索正链，通过-m/--max-mismatch 来允许错配 -d, --degenerate 包含简并碱基的pattern和motif（简并碱基：一个符号代替某两个或者更多碱基，如编译丙氨酸的可以有4个密码子：GCU\GCC\GCA\GCG,这时生物学上为了方便，用字母N指代UCAG四个碱基，故说编译丙氨酸的密码子是GCN，其中N就是简并碱基。） -i, --ignore-case 忽视大小写 -v, --invert-match 输出不匹配此模式的内容 -m, --max-mismatch int 通过序列匹配时允许的最大错配 -p, --pattern strings 匹配模式，支持连续写多个模式，匹配任一模式即输出。 -f, --pattern-file string 支持匹配模式写到一个文件中，如要提取的序列ID。 -R, --region string 搜索特定区域序列，e.g 1:12 for first 12 bases, -12:-1 for last 12 bases -r, --use-regexp 使用正则表达式，必须加入此参数,如^匹配首端。同-p联合使用。 举例： seqkit grep -s -r -i -p ^atg cds.fa＃选取正链中有起始密码子的序列 seqkit grep -f ID.txt test.fa > new.fa＃根据ID提取序列 seqkit grep -s -d -i -p TTSAA＃简并碱基使用。S 代表C or G. seqkit grep -s -R 1:30 -i -r -p GCTGG＃＃匹配限定到某区域 usage:seqkit replace(通过正则表达式来代替名字或序列) -s, --by-seq 代替seq -i, --ignore-case 忽视大小写 -K, --keep-key -p, --pattern string 搜索正则表达式 -r, --replacement string 替换物从两个配对端读数的文件提取配对的reads ## 首先创造两个不平衡的PE reads文件，整个过程不报错，但是不懂具体含义？？？ $ seqkit rmdup duplicated-reads.fq.gz \ | seqkit replace --pattern " .+" --replacement " 1" \ | seqkit sample --proportion 0.9 --rand-seed 1 --out-file read_1.fq.gz $ seqkit rmdup duplicated-reads.fq.gz \ | seqkit replace --pattern " .+" --replacement " 2" \ | seqkit sample --proportion 0.9 --rand-seed 2 --out-file read_2.fq.gz ## overview $ seqkit stat read_1.fq.gz read_2.fq.gz file format type num_seqs sum_len min_len avg_len max_len read_1.fq.gz FASTQ DNA 9,033 912,333 101 101 101 read_2.fq.gz FASTQ DNA 8,965 905,465 101 101 101 ## sequence headers $ seqkit head -n 3 read_1.fq.gz | seqkit seq --name SRR1972739.1 1 SRR1972739.3 1 SRR1972739.4 1 $ seqkit head -n 3 read_2.fq.gz | seqkit seq --name SRR1972739.1 2 SRR1972739.2 2 SRR1972739.3 2 ## 首先提取两个文件序列的ID，并计算它们的交集 $ seqkit seq --name --only-id read_1.fq.gz read_2.fq.gz \ | sort | uniq -d > id.txt $ # number of IDs wc -l id.txt 8090 id.txt ## 然后用id.txt来提取reads $ seqkit grep --pattern-file id.txt read_1.fq.gz -o read_1.f.fq.gz $ seqkit grep --pattern-file id.txt read_2.fq.gz -o read_2.f.fq.gz ## 用md5sum来检查两个文件的IDs是否是相同的 $ seqkit seq --name --only-id read_1.f.fq.gz > read_1.f.fq.gz.id.txt $ seqkit seq --name --only-id read_2.f.fq.gz > read_2.f.fq.gz.id.txt $ md5sum read_*.f.fq.gz.id.txt 537c57cfdc3923bb94a3dc31a0c3b02a read_1.f.fq.gz.id.txt 537c57cfdc3923bb94a3dc31a0c3b02a read_2.f.fq.gz.id.txt ##sort 一下 $ gzip -d -c read_1.f.fq.gz \ | seqkit fx2tab \ | sort -k1,1 -T . \ | seqkit tab2fx \ | gzip -c > read_1.f.sorted.fq.gz $ gzip -d -c read_2.f.fq.gz \ | seqkit fx2tab \ | sort -k1,1 -T . \ | seqkit tab2fx \ | gzip -c > read_2.f.sorted.fq.gz将两个fasta文件连接成一个 $ cat 1.fa >seq1 aaaaa >seq2 ccccc >seq3 ggggg $ cat 2.fa >seq3 TTTTT >seq2 GGGGG >seq1 CCCCC 一行命令 $ seqkit concat 1.fa 2.fa >seq1 aaaaaCCCCC >seq2 cccccGGGGG >seq3 gggggTTTTT |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |