9. seqtk seqkit gtftk 总结

1. 背景

  在前面小节我们使用了这些软件,因为混合使用比较让人混乱,这里总结理清楚一下.

2. seqtk

  功能总览如下图所示.
image

2.1 seq

  这个功能主要是对\(.fasta\)和\(.fastq\)格式的文件进行格式化.
image

  • \(-l\)
      主要是让序列每行显示多少个碱基
#每行显示60个氨基酸
seqtk seq -l 60 atha.fasta | less -SN
  • -L
      将短于\(x\)的序列删掉.
# 将短于100的氨基酸序列删除
seqtk seq -L 100 atha.fasta | less -SN

image

  • -U
      将序列里的小写字母转化为大写,这个应该不用演示了.

2.2 Sample

  当我们序列条数过多,全部用于实验会让代码速度减慢,所以有时我们会将序列抽样
image
  从上图理解,格式是seqtk sample in.fa 分数|具体数字,也就是说我们可以抽取多少含量,也可以抽具体数字.

seqtk sample atha.fasta 2

image
  如果我们不改变随机种子,那么每次随机都会是一样的结果.

  • -s
      改变随机种子

2.3 subseq

  用此指令提取序列.
image
  可以观察到第一个参数是源文件,第二个参数是对应键名文件,我们根据\(name.list\)去提取文件.

seqtk subseq genome.fa name.list | less -N

我们可以改变name.list的文件内容,让\(subseq\)提取不同位置的碱基.代码保持不变,获得的碱基不同了.
image
image
  这里有一个点需要科普.我们在文件里列的是1~10,但最后展示出来的是2~10.虽然我们起的名是\(name.list\)但是实际标准输入文件是bed格式的文件.在bed文件中,染色体序列的序号是从\(0\)开始,而输出文件是从\(1\)开始,这就造成了两者偏差.

3. seqkit

  \(seqkit\)的很多功能与\(seqtk\)重叠.但是还是有一些功能用\(seqkit\)比较多.

3.1 faidx

  此指令用于构建Index.创建 FASTA 索引文件并提取子序列.这个指令类似于samtools faidx,网上关于这个指令的资料没多少,但是samtools的资料蛮多.下面是参考\(samtools\)的资料.
  samtools faidx(seqkit faidx也可以) 能够对fasta 序列建立一个后缀为.fai 的文件，根据这个.fai 文件和原始的fastsa文件， 能够快速的提取任意区域的序列.该命令对输入的fasta序列有一定要求：对于每条序列，除了最后一行外， 其他行的长度必须相同.生成的\(fai\)文件共5列,\t分隔.写着写着感觉我好像在上一节写过,算了不管了,接下来搞懂makeblastdb.
image

3.2 shuffle

  显然是用于打乱序列.默认情况下,所有记录都会被读入内存.对于 FASTA 格式，使用标志 -2 (--two-pass) 来减少内存使用.顺带一提,不支持FASTQ的格式.
image

seqkit shuffle atha.fasta | less -SN

上面都是fasta的文件处理,下面是gff、gtf的文件处理
上面都是fasta的文件处理,下面是gff、gtf的文件处理
上面都是fasta的文件处理,下面是gff、gtf的文件处理

4. gffread

  这个软件在现时期已经变成了\(gtf\)和\(gff\)文件转换器.虽然它还有其他功能.

4.1 gff2gtf

# GFF转换GTF
gffread input.gff3 -T -o out.gtf
#GTF转换GFF3
gffread input.gtf -o out.gff3

参数列表:

-g   序列文件，即GFF/GTF文件第一列ID对应的序列文件。
-i   丢弃掉内含子大于的转录本（mRNA/transcript）
-r   起始和终止位置，填写示例100.10000即为输出与100到10000有重叠的所有转录组，也可以限制序列ID及链，填写示例：+Chr1:100..10000。
-R   丢弃掉此范围的转录本，与-r相反。
-U   丢弃掉 single-exon的转录本
-C   丢低调无CDS的转录本。
-V   丢弃掉含有移码突变的转录本。
-H   如果使用了-V，则重新检查并调整内含子相位，避免由于翻译起始位点选择的位置不对导致移码突变的产生。
-B   如果使用了-V， 对于单外显子基因，则重新检查相反的链，是否存在移码突变。
-N   丢弃掉多外显子基因剪接位点不是常见的 GT-AG, GC-AG or AT-AC序列。
-J   丢弃掉没有起始密码子或者终止密码子的转录本，仅保留有完整编码框的转录本。
--no-pseudo:  过滤掉含有 'pseudo' 的注释信息
-M/--merge :  合并完全相同的或者存在包含关系的转录本。
-d:使用 -M  将合并信息输出到文件中
--cluster-only: 类似于 --merge 但是不合并转录本
-K   对于-M 选项：also collapse shorter, fully contained transcripts
      with fewer introns than the container
-Q   对于-M 选项：移除包含关系的转录本的限制条件：多外显子转录本将会合并，如果他们内含子位置完全一样，单外显子转录本只需要有80%一样即可合并。
--force-exons:  使GFF features的最小层级为exon
-E   对于重复的 ID或者 GFF/GTF 其他潜在的格式问题给出警告信息。
-Z   将内含子小于4 bp的邻近的两个外显子合并为一个。
-w   输出每个转录本的外显子序列
-x   输出CDS序列
-W   对于 -w 和 -x 选项，输出外显子位置坐标到FASTA序列的ID中
-y   输出蛋白质序列
-L   将Ensembl GTF 转换为 GFF3 conversion (implies -F; should be used with -m)
-o   输出"filtered" 后的GFF文件 。
-T -o  参数将输出 GTF格式。

5. pygtftk

  为什么网上没有这个包的资料...总之翻译了一下官网的话.Python GTF 工具包 (pygtftk) 包旨在简化 GTF/GFF2.0 文件（基因传输格式）的处理。目前不支持 GFF3 文件格式。 pygtftk 包与 Python >=3.5、<3.7 兼容，并依赖于 libgtftk(一个用 C 编写的函数库).这里注意一下,\(gtf\)必须是\(ensembl\)格式的(\(gtf\)的\(ensemble\)格式是在\(feature\)列有\(gene\)、\(transcript\)列,如果没有就是常规的).
  本人\(python\ 3.9\)不知道为啥也能用.

5.1 conversion

  • convert_ensembl
      将 GTF 文件转换为 ensembl 格式。本质上是添加“转录本”/“基因”功能.懒了直接复制的上一小节.
gtftk convert_ensembl -i temp.gtf > SWO_genes.gtf
#-i, --inputfile              Path to the poor BED file to would like to behave as if it was a GTF.(default: <stdin>)
#-o, --outputfile             Output file. (default: <stdout>)
  • tabulate
      从 GTF 中提取键/值并将其转换为表格格式。当请求坐标时，它们将以基于 1 的格式提供.
gtftk nb_transcripts -i SWO_genes.gtf | gtftk select_by_key -k feature -v gene | gtftk tabulate -k gene_id,nb_tx -s "|" | less -SN

image

5.2 information

  这部分的指令主要用于提取信息.

  • count
      计算特征的数量(转录本、基因、外显子、内含子).
Usage: gtftk count [-i GTF] [-o TXT] [-d header] [-t TEXT] [-h] [-V [verbosity]] [-D] [-C] [-K tmp_dir] [-A] [-L logger_file] [-W write_message_to_file]

$ gtftk count -i simple.gtf -t example_gtf
gene	10	example_gtf
transcript	15	example_gtf
exon	25	example_gtf
CDS	20	example_gtf
  • nb_transcripts
      计算每个基因的转录本数量.我就不放参数了,直接去看官方文档,详细多了.
gtftk nb_transcripts -i SWO_genes.gtf | less -SN

得到的输出如下所示,我们可以发现每个\(feature\)为\(gene\)的行末都有一个\(nb_tx\),这是记录这个基因有多少个转录本.
image

  • nb_exons
      计算外显子的数量并将其添加为新的键/值。如果需要，输出也可以是文本格式.
$ gtftk nb_exons -i simple.gtf | head -n 5
chr1	gtftk	gene	125	138	.	+	.	gene_id "G0001";
chr1	gtftk	transcript	125	138	.	+	.	gene_id "G0001"; transcript_id "G0001T002"; nb_exons "1";
chr1	gtftk	exon	125	138	.	+	.	gene_id "G0001"; transcript_id "G0001T002"; exon_id "G0001T002E001";
chr1	gtftk	CDS	125	130	.	+	.	gene_id "G0001"; transcript_id "G0001T002"; ccds_id "CDS_G0001T002";
chr1	gtftk	transcript	125	138	.	+	.	gene_id "G0001"; transcript_id "G0001T001"; nb_exons "1";

image

  • feature_size
      获取 GTF 中包含的\(feature\)的大小和限制（开始/结束）。如果请求bed格式，则返回bed格式的限制和大小作为分数。否则输出 GTF 文件，其中"feat_size"作为新键，大小作为值.我人话翻译一下,就是计算这个\(feature\)的序列大小,默认计算\(feature=transcript\)的碱基个数.
    image

5.3 select

  从源文件中选择某些基因.

  • short_long
      获取每个基因的最短或最长转录本.
    参数：
  1. -i, --inputfile GTF 文件的路径。 默认为 STDIN(默认值：stdin).
  2. -o, --outputfile 输出文件。 （默认值：<标准输出>）
  3. -l, --longs 取每个基因最长的转录本（默认值：False）
  4. -g, --keep-gene-lines 将基因行添加到输出（默认值：False）

  这个输出文件有点长,不适合放在博客上展示.

  • select_by_key
      根据键和值从 gtf 中提取行
gtftk nb_transcripts -i SWO_genes.gtf | gtftk select_by_key -k feature -v gene | less -SN
#选取 key=feature value=gene的行

image

5.4 sequence

  获取序列,从下文看,是根据下载的基因组文件和基因组注释文件,将对应部件的序列列出来,同时也可以提供序列开头和末尾的信息.

  • get_tx_seq
      获取fasta格式的\(transcript\)序列
gtftk get_tx_seq -g genome.fa  -i SWO_genes.gtf less -SN
#  -g, --genome fasta 格式的基因组。接受带有通配符的路径
#  -i, --inputfile GTF 文件的路径。默认为 STDIN（默认值：<stdin>）
#  -l, --label 标题的一组键.(默认：feature、transcript_id、gene_id、seqid、start、end)

image

  • get_feat_seq
      获取特征序列(例如外显子、UTR…).这里直接复制了
$ gtftk get_feat_seq -i simple.gtf -g simple.fa  -l feature,transcript_id,start -t  exon -n | head -10
index file simple.fa.fai not found, generating...
>exon|G0001T002|124
cccccgttacgtag
>exon|G0001T001|124
cccccgttacgtag
>exon|G0002T001|179
ggccttatta
>exon|G0003T001|49
caagc
>exon|G0003T001|56
taatt
posted @ 2023-09-29 14:43  acmloser  阅读(680)  评论(0)  编辑  收藏  举报