搜索
查看: 2654|回复: 0

GATK之FastaStats和CountReads

[复制链接]

18

主题

56

帖子

403

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
403
发表于 2017-5-4 20:25:27 | 显示全部楼层 |阅读模式
本帖最后由 hoptop 于 2017-5-4 20:30 编辑

GATK的本职工作是Variant calling,但是就像我之前所说的,它作为基因组分析工具箱,还是有很多其他工具,今天学习的是诊断和质量控制工具的其中两个:CountReads,FastaStats。



FastaStats
功能: 计算参考基因组的基本统计值
分类: 诊断和质量控制工具
概要: 主要就是统计碱基总数,和常规的碱基数(ATCG)

吐槽:功能还真是简单。。

输入: FASTA参考文件
输出: 结果要么输出到屏幕,要么是输出到(-o)到文件中。
使用案例:我看了一下拟南芥的参考基因组
[Bash shell] 纯文本查看 复制代码
java -jar ~/biosoft/GenomeAnalysisTK.jar -T FastaStats \
    -R TAIR10.fa
输出:
Total bases   119667750Regular bases 119481543






CountReads
功能: 计算reads数
分类: 诊断和质量控制工具
概要: 最好和--read-filter合用,这样可以了解下符合特定标准的reads数
输入: 一个或多个BAM文件
输出: 结果会输出到屏幕(标准输出)上,毕竟是用来确定阈值的,也不需要一定要输出到文件中

案例:不加参数--read-filter

[AppleScript] 纯文本查看 复制代码
java -jar ~/biosoft/GenomeAnalysisTK.jar -T CountReads \
    -R $work/database/TAIR10/TAIR10.fa  -I BC_bg_reads.sorted.bam
输出:
CountReads - CountReads counted 55080781 reads in the traversal


--read-filter/-rf后面可以接很多的选项,官方文档列出了如下内容:

NameSummary
BadCigarFilterFilter out reads with wonky CIGAR strings
BadMateFilterFilter out reads whose mate maps to a different contig
CountingFilteringIterator.CountingReadFilter
DuplicateReadFilterFilter out duplicate reads
FailsVendorQualityCheckFilterFilter out reads that fail the vendor quality check
HCMappingQualityFilterFilter out reads with low mapping qualities for HaplotypeCaller
LibraryReadFilterOnly use reads from the specified library
MalformedReadFilterFilter out malformed reads
MappingQualityFilterFilter out reads with low mapping qualities
MappingQualityUnavailableFilterFilter out reads with no mapping quality information
MappingQualityZeroFilterFilter out reads with mapping quality zero
MateSameStrandFilterFilter out reads with bad pairing (and related) properties
MaxInsertSizeFilterFilter out reads that exceed a given insert size
MissingReadGroupFilterFilter out reads without read group information
NoOriginalQualityScoresFilterFilter out reads that do not have an original quality quality score (OQ) tag
NotPrimaryAlignmentFilterFilter out read records that are secondary alignments
OverclippedReadFilterFilter out reads that are over-soft-clipped
Platform454FilterFilter out reads produced by 454 technology
PlatformFilterFilter out reads that were generated by a specific sequencing platform
PlatformUnitFilterFilter out reads with blacklisted platform unit tags
ReadGroupBlackListFilterFilter out reads matching a read group tag value
ReadLengthFilterFilter out reads based on length
ReadNameFilterOnly use reads with this read name
ReadStrandFilterFilter out reads based on strand orientation
ReassignMappingQualityFilterSet the mapping quality of all reads to a given value
ReassignOneMappingQualityFilterSet the mapping quality of reads with a given value to another given value
ReassignOriginalMQAfterIndelRealignmentFilterRevert the MQ of reads that were modified by IndelRealigner
SampleFilterOnly use reads belonging to a specific sample
SingleReadGroupFilterOnly use reads from the specified read group
UnmappedReadFilterFilter out unmapped reads

使用的时候记住,对于XXXXFilter,你需要点每个过滤选项进去看下具体用法,比如说我需要过滤比对质量低于20的,我发现MappingQualityFilter的用法解释是
[AppleScript] 纯文本查看 复制代码
java -jar GenomeAnalysisTk.jar \
         -T HaplotypeCaller \
         -R reference.fasta \
         -I input.bam \
         -o output.vcf \
         -rf MappingQuality \
         -mmq 15


那么我前面的代码就可以改为

[AppleScript] 纯文本查看 复制代码
java -jar ~/biosoft/GenomeAnalysisTK.jar -T CountReads -R $work/database/TAIR10/TAIR10.fa  -I BC_bg_reads.sorted.bam 
\-rf MappingQuality -mmq 15
结果:
8190205 reads were filtered out during the traversal out of approximately 55080781 total reads (14.87%)

果然加上比对质量后,就有一些要被过滤掉呢。
基本上GATK每一个工具都有一些默认过滤选项,你可以在每个工具的额外信息部分进行了解。



今天就介绍这两个工具,顺便提了GATK的过滤功能。白了个白,每天学个GATK的任务完成。

PS: 欢迎关注公众号 生信媛, 里面有我写的很多文章呢。







上一篇:5种场合需要用action button
下一篇:为什么不能用t检验取代方差分析
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-14 22:13 , Processed in 0.034304 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.