基因组变异和 NGS 读取的命令行操作
项目描述
瓦伦斯
一组用于处理基因组变异和下一代测序读取的 Python 工具。对于大型数据集不是特别快。重点是将您需要的 BAM 和 VCF 提取到 CSV 文件中以供进一步分析。
- varlens 变体
组合、注释和过滤来自 VCF 或 CSV 文件的变体。可用的注释包括基因、变异效应、周围序列上下文、来自指定 BAM 文件的支持读取计数以及突变肽的 MHC I 结合亲和力预测。
- varlens 读取
显示、过滤和复制来自 SAM/BAM 文件的读数。samtools 视图的部分替换。
- varlens等位基因支持
计数支持 BAM 文件中指定位点的每个等位基因的读数。
安装
从PyPI安装：
pip install varlens
或者来自 git checkout：
pip install .
要运行测试：
nosetests .
要构建文档（仅此自述文件和命令行工具帮助）：
pip install -e . pip install Sphinx cd docs make clean setup rst html
文档将被写入_build/html目录。
varlens 变体
给定来自一个或多个 VCF 或 CSV 文件的变体，应用过滤器，添加其他列，然后输出到 CSV。
目前我们只能输出到 CSV，不能输出到 VCF。
通过指定“–include-XXX”形式的选项，例如“–include-gene”，可以为每个变体添加许多有用的注释。查看详细帮助（使用 -h 运行）。
例子
打印在两个 VCF 文件中找到的变体的基本信息。请注意，在这两个文件中找到的变体都列在一行中，并且“来源”列列出了在其中找到每个变体的文件：
$ varlens-variants test/data/CELSR1/vcfs/vcf_1.vcf test/data/CELSR1/vcfs/vcf_2.vcf genome,contig,interbase_start,interbase_end,ref,alt,sources GRCh37,22,21829554,21829555,T,G,1.vcf GRCh37,22,46931059,46931060,A,C,1.vcf GRCh37,22,46931061,46931062,G,A,1.vcf 2.vcf GRCh37,22,50636217,50636218,A,C,1.vcf GRCh37,22,50875932,50875933,A,C,1.vcf GRCh37,22,45309892,45309893,T,G,2.vcf
与上述相同，但包括提供 varcode 变体效果注释的附加列以及变体重叠的基因，并写入文件：
$ varlens-variants test/data/CELSR1/vcfs/vcf_1.vcf test/data/CELSR1/vcfs/vcf_2.vcf \
--include-effect \
--include-gene \
--out /tmp/result.csv
Wrote: /tmp/result.csv
$ cat /tmp/result.csv
genome,contig,interbase_start,interbase_end,ref,alt,sources,effect,gene
GRCh37,22,21829554,21829555,T,G,1.vcf,non-coding-transcript,PI4KAP2
GRCh37,22,46931059,46931060,A,C,1.vcf,p.S670A,CELSR1
GRCh37,22,46931061,46931062,G,A,1.vcf 2.vcf,p.S669F,CELSR1
GRCh37,22,50636217,50636218,A,C,1.vcf,intronic,TRABD
GRCh37,22,50875932,50875933,A,C,1.vcf,splice-acceptor,PPP6R2
GRCh37,22,45309892,45309893,T,G,2.vcf,p.T214P,PHF21B
打印支持来自指定 BAM 的参考/变体/其他等位基因的读数的计数，仅计算映射质量 >= 10 的读数：
$ varlens-variants test/data/CELSR1/vcfs/vcf_1.vcf \
--include-read-evidence \
--reads test/data/CELSR1/bams/bam_1.bam \
--min-mapping-quality 10
genome,contig,interbase_start,interbase_end,ref,alt,sources,num_alt,num_ref,total_depth
GRCh37,22,21829554,21829555,T,G,vcf_1.vcf,0,0,0
GRCh37,22,46931059,46931060,A,C,vcf_1.vcf,0,216,320
GRCh37,22,46931061,46931062,G,A,vcf_1.vcf,0,321,321
GRCh37,22,50636217,50636218,A,C,vcf_1.vcf,0,0,0
GRCh37,22,50875932,50875933,A,C,vcf_1.vcf,0,0,0
varlens 读取
过滤来自一个或多个 BAM 的读取并输出 CSV 或新的 BAM。
可以指定基因座和 VCF 文件，在这种情况下，读数会被过滤以与指定的基因座或变体重叠。
例子
打印 BAM 中读取的基本字段：
$ varlens-reads test/data/CELSR1/bams/bam_0.bam query_name,reference_start,reference_end,cigarstring HISEQ:142:C5822ANXX:3:2116:16538:101199,46929962,46930062,100M HISEQ:142:C5822ANXX:3:1106:18985:32932,46929964,46930064,100M HISEQ:142:C5822ANXX:3:2201:21091:67220,46929966,46930066,100M HISEQ:142:C5822ANXX:4:1304:5363:12786,46929966,46930066,100M HISEQ:142:C5822ANXX:4:1104:9008:85114,46929969,46930069,100M HISEQ:142:C5822ANXX:3:2304:9921:94828,46929970,46930070,100M HISEQ:142:C5822ANXX:3:2211:6266:74633,46929973,46930073,100M HISEQ:142:C5822ANXX:3:1305:8982:42729,46929974,46930074,100M HISEQ:142:C5822ANXX:4:2316:5630:7371,46929978,46930078,100M ...
与上面相同，但仅过滤在 (-) 链上对齐的读取，写入文件而不是标准输出，并且还在输出中包含映射质量和测序碱基：
$ varlens-reads test/data/CELSR1/bams/bam_0.bam \
--is-reverse \
--field mapping_quality query_alignment_sequence \
--out /tmp/result.csv
Wrote: /tmp/result.csv
$ head /tmp/result.csv
query_name,reference_start,reference_end,cigarstring,mapping_quality,query_alignment_sequence
HISEQ:142:C5822ANXX:3:2116:16538:101199,46929962,46930062,100M,60,CATGATCTGGGCATTAGGGCCTTCATCAGGGTCGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCG
HISEQ:142:C5822ANXX:3:1106:18985:32932,46929964,46930064,100M,60,TGATCTGGGCATTAGGGCCTTCATCAGGGTCGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTC
HISEQ:142:C5822ANXX:4:1104:9008:85114,46929969,46930069,100M,60,TGGGCATTAGGGCCTTCATCAGGGTCGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCT
HISEQ:142:C5822ANXX:4:1202:18451:91174,46929979,46930079,100M,60,GGCCTTCATCAGGGTCGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCTCAAACATGGG
HISEQ:142:C5822ANXX:3:1211:18522:54773,46929987,46930087,100M,60,TCAGGGTCGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCTCAAACATGGGGGCATTGT
HISEQ:142:C5822ANXX:3:2114:19455:45093,46929987,46930087,100M,60,TCAGGGTCGTTAGCACGAATCTTTGCCACCGCCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCTCAAACATGGGGGCATTGT
HISEQ:142:C5822ANXX:4:2115:9153:21593,46929994,46930094,100M,60,CGTTAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCTCAAACATGGGGGCATTGTCATTAAT
HISEQ:142:C5822ANXX:4:1212:15644:87227,46929995,46930095,100M,60,GTTAGCACGTATGTTTGCCACCACCGACCCCACTGAGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTGCTTCTCAAACATGGGGGCAGTGTCATTAATG
HISEQ:142:C5822ANXX:3:1103:4717:26369,46929997,46930097,100M,60,TAGCACGAATCTTTGCCACCACCGACCCCACTGGGTTGTTCTCCTCAACAAACAGCTCCAGTTCGTCCTTCTCAAACATGGGGGCATTGTCATTAATGTC
编写一个 bam 文件，其中包含映射质量 >=30 且与某个位点重叠的读取：
$ varlens-reads test/data/CELSR1/bams/bam_0.bam \
--min-mapping-quality 30 \
--locus 22:46932040-46932050 \
--out /tmp/result.bam
编写一个 bam 文件，其中包含来自 VCF 的读取重叠变体：
$ varlens-reads test/data/CELSR1/bams/bam_0.bam \
--variants test/data/CELSR1/vcfs/vcf_1.vcf \
--out /tmp/result.bam
仅以 csv 格式打印 BAM 的标题：
$ varlens-reads test/data/CELSR1/bams/bam_0.bam --header
varlens等位基因支持
给定一个或多个 BAM 和一些要考虑的基因组位点，编写一个 csv 文件，给出每个 BAM 在每个位点支持每个等位基因的读取计数。
要考虑的基因组位点可以通过基因座（–locus 选项）或通过一个或多个 VCF 文件指定。
该命令输出的位置是在基间坐标中，即从 0 开始，包括第一个索引，不包括第二个索引（与 VCF 文件中使用的基于一的包含坐标相反）。
例子
varlens-allele-support \
--reads test/data/CELSR1/bams/bam_1.bam \
--locus 22:46931061 22:46931063
source,contig,interbase_start,interbase_end,allele,count
bam_1.bam,22,46931060,46931061,,1
bam_1.bam,22,46931060,46931061,G,329
bam_1.bam,22,46931062,46931063,A,327
bam_1.bam,22,46931062,46931063,AC,1
bam_1.bam,22,46931062,46931063,AG,2
坐标系注意事项
varlens 在内部使用基于 0 的半开坐标。许多工具（包括 samtools 和 VCF 文件）使用包含 1 的坐标。当我们使用基于 0 的半开坐标时，我们使用术语“interbase”来尽量减少混淆，而当我们使用基于 1 的包容性坐标时，我们会使用术语“inclusive”。
一个特别棘手的地方是在命令行上使用例如--locus chr22:43243-43244指定基因座时。为了与最常见的其他工具保持一致，当您指定像 chr22:10-20这样的轨迹时，我们将其解释为基于 1 的包含坐标。要指定从 0 开始的半开坐标，请使用以下语法：chr22/11-20（即斜杠而不是冒号）。
有关坐标系的更多详细信息，请参阅此博客文章 。
