samtools、bedtools、fastxtoolkit.ppt

samtools、bedtools、fastx-toolkit,北京大学生命科学学院叶永鑫2011年5月23日星期一,outline,samtools操作sam/bam文件(mapping)的工具附加bcftools操作vcf/bcf文件(snp/indelcalling)的工具bedtools操作bed文件(blockfeatureannotation)的工具fastx-toolkit操作fasta/fastq文件的工具,samtools,sam文件格式,sequencealignment/mapformat存储mapping的结果的tab分隔的文本文件headersection(optional)开头，例如：sqsn:refln:45alignmentsection11个必需的字段，例如：r001163ref7308m2i4m1d3m=3739ttagataaaggatactg*,headersection,hdheaderlinevn*formatversionsosortingorderofalignmentssqreferenceseq.dictionarysn*ln*ref.seq.lengthasgenomeassemblyidentifierm5md5checksumspspeciesururi,headersection,rgreadgroupid*readgroupidentifiercnseq.centerdsdescriptiondtdateoftherunfofloworderkskeysequencelblibrarypgprogramspipredictedmedianinsertsizeplplatform/technologypuplatformunitsmsample,headersection,pgprogramid*programrecordidentifierpnprogramnameclcommandlinepppreviouspg-idvnprogramversioncocomment(one-linetext),alignmentsection-mandatoryfields,1qnamequerytemplatename2flagbitwiseflag3rnameref.seq.name4pos1-basedleftmoseposition5mapqmappingquality6cigarcigarstringmmatch/mismatchiinsertionddeletion（相对reference）ssoftclippinghhardclippingppaddingnskip=matchxmismatch,alignmentsection-mandatoryfields,7mrnmmateref.sequencename(=ifsameasrname)8mpos1-basedmateposition9isizeinferredinsertsize10seqquerysequence11qualqueryquality(asciiofphred-scaledbasequality),alignmentsection-mandatoryfields-flag,0 x001pthereadispairedinsequencing0 x002pthereadismappedinaproperpair0 x004uthequeryseq.itselfisunmapped0 x008uthemateisunmapped0 x010rstrandofthequery(1forreverse)0 x020rstrandofthemate0 x0401thereadisthefirstreadinapair0 x0802thereadisthesecondreadinapair0 x100sthealignmentisnotprimary0 x200fthereadfailsplatformqualitychecks0 x400dthereadiseitherapcroropticaldup.,alignmentsection-optionalfields,tag:type:valuecscolorreadseq.cqcolorreadqualitymdstringformismatchingpositionaimtoachievesnp/indelcallingwithoutlookingatthereferencenmeditdistancetothereferncergreadgroupmatchrg-idifrgexists,bam文件格式,binarysamcompressedinthebgzfformat查看bam文件内容：samtoolsview-haln.bam|lesssbam可被排序：samtoolssortaln.bamaln.sorted排序的bam可建索引：samtoolsindexaln.sorted.bam,samtools,sam、bam转换view、import查看sam、bam文件view、tview排序bam文件sort建索引faidx、index、tabix*,samtools,初步统计flagstats、idxstats、depth*改动bam文件reheader、merge、cat*校正bam文件fixmate、calmd、rmdupcallsnp/smallindelpileup、mpileup,samtools,samtools设计为可以在流（管道）上执行，在需要file处填-就可表示stdin或stdout。命令行未跟足够参数时就会显示samtools各程序的简要帮助。eg.samtoolssamtoolsviewsamtoolssort,sam、bam转换,sam-bamsamtoolsimportref.fa.faiin.samout.bamifin.samw/sq:samtoolsview-bsin.samout.bamifin.samw/osq:samtoolsfaidxref.fasamtoolsview-btref.fa.faiin.samout.bamorsamtoolsview-btref.fain.samout.bambam-samsamtoolsview-hin.bamout.sam,samtoolsview,samtoolsviewoptionsin.bam/samregion1.-?longerhelp-boutputinbamformat-cinstead,onlycount-fintonlyoutputalignmentswithallbitsinflag0-fintskipalignmentswithflagbits0-hincludetheheaderinoutput-houtputtheheaderonly-lstronlyoutputreadsinlibrarystrnull,samtoolsview,-ofileoutputfilestdout-qintskipalignmentswithmapqualitysmallerthanint0-rstronlyoutputreadsinreadgroupstrnull-rfileoutputreadsinreadgroupslistedinfilenull-sinputisinsam,ifsqisabsent,-tor-tisrequired-tfile指定referencesequencefile.若尚未建立索引，会自动先建立索引，类似ref.fa.fai,samtoolsview,-tfiletab-delimitedfile,andlength;canbetheresultantindexfileref.fa.faibysamtoolsfaidxref.fa-uoutputuncompressedbam,forpipe-xoutputflaginhex(eg.0 x00)-xoutputflaginstringregion(1-based)chr2chr2:1000(startfrom1000)chr2:1,000-2,000(include2000),samtoolsview,eg.samtoolsview-bsin.samout.bamsamtoolsview-btref.fa.faiin.samout.bamsamtoolsview-btref.fain.samout.bamsamtoolsview-hin.bamout.samsamtoolsviewex1.bamseq1:1-10|less-snote:randomretrievalneedindexedsamtoolsview-xex1.bam|less-ssamtoolsview-xex1.bam|less-s,对bam文件排序和建索引,samtoolssort可以对bam文件进行排序。samtoolsindex可以对已排序的bam文件建立索引，从而使得对该bam文件的randomretrieval成为可能。eg.samtoolssortaln.bamaln.sorted（产生已排序的aln.sorted.bam文件）samtoolsindexaln.sorted.bam（产生aln.sorted.bam.bai索引文件）,samtoolssort,samtoolssortoptionsin.bamout.prefix-ooutputthefinalalignmenttothestandardoutput即把排序的bam文件输出到stdout；仍然需要在命令行指定out.prefix，否则不能执行，但不产生out.prefix.bam-nsortbyreadnamesratherthanbychromosomalcoordinates-mintapproximatelythemaximumrequiredmemory,建立索引,建立索引的目的是使randomretrieval成为可能samtoolsindexaln.bam为已排序的aln.bam建索引，将创建aln.bam.baisamtoolsfaidxref.fastaregion1.为referencefasta建索引，将创建ref.fasta.fai,samtoolstview,samtoolstviewin.sorted.bamref.fastain.sorted.bam需要先用samtoolsindex建索引ref.fasta可不用先建索引；若无索引，samtoolstview会自动去建立ref.fasta.fai索引insamtoolstview:?helpggotoaregionchrm:1000;seq1:10=1000(ifsamereferencesequence)qexit,samtoolstview,h,j,k,lindexfilefile.faiwillbecreatedifabsent-gfloatpriorofanindelbetweentwohaplotypes(for-c)0.00015-ionlyshowlines/consensuswithindels-iintphredprob.ofanindelinsequencing/prep.40-lfilelistofsitesatwhichpileupisoutput,samtoolspileup,-mintfilterreadswithflagcontainingbitsinint1796(0 x704,usfd)-mintcapmappingqualityatint60-nint#haplotypesinthesample(for-c)(=2)2-qintminbasequality13-rfloatpriorofadifferencebetweentwohaplotypes(for-c)0.001-stheinputisinsam,samtoolspileup,-tfilelistofreferencenamesandsequencelengths(force-s)-tfloatthethetaparameter(errordependencycoefficient)inmaqconsensuscallingmodel(for-c)0.83-vprintvariantsonly(for-c)note:pileupisdeprecated,pleaseusempileup,samtoolspileup,defaultoutput:pileupformat每一行都是一个genomicposition，包含以下6列：1chromosomename2coordinate3referencebase4readcoverage/depth5readbases6alignmentmappingqualities,samtoolspileup,readbases.foramatchtotheref.onthe+strand,foramatchtotheref.onthe-strandflt.txtpleaseremembertosetthe-dtosetamaximumreaddepthaccordingtotheaveragereaddepth(2).inwholegenomeshotgun(wgs)resequencing,snpswithexcessivelyhighreaddepthareusuallycausedbystructuralvariationsoralignmentartifactsandshouldnotbetrusted.awk($3=*0forref,1foralt1,.plphred-scaledgenotypelikelihoodorder:00,01,11,02,12,22,.gqconditionalgenotypephredquality10+sample(s)individualgenotypeinformationdefinedbyformat,bcftools,bcftoolsview查看，转换文件格式index为bcf建索引cat连接bcfldcomputeall-pairr2ldpaircomputer2betweenrequestedpairs,bcftoolsview,bcftoolsviewoptionsin.bcfregion-boutputbcfinsteadofvcf-csnpcalling(force-e)-elikelihoodbasedanalyses-gsuppressallindividualgenotypeinfo.-gcallgenotypesatvariantsites(force-c)-iskipindels-lfilelistofsites(chrpos)orregions(bed)tooutputall-pstrtypeofprior:full,cond2,flatfull-voutputpotentialvariantsitesonly(force-c),vcfutils.pl,vcfutils.plvarfiltervcf2fqsubsamlistsamfillacqstatshapmap2vcfucscsnp2vcf,vcfutils.plvarfilter,vcfutils.plvarfilteroptionsin.vcf-aintminnumberofalternatebases2-dintmaxreaddepth10000000-dintminreaddepth2-pprintfilteredvariantstostderr-qintminrmsmappingqualityforsnps10-wintwindowsizeforfilteringadjacentgaps10-wintsnpwithintbparoundagaptobefiltered10,vcfutils.plvcf2fq,vcfutils.plvcf2fqoptionsall-site.vcf输入samtoolsmpileup-uf|bcftoolsview-cg的结果，转换为fastq格式-dintmindepth3-dintmaxdepth100000-lintindelfilterwindowsize5-qintminrmsmappingquality10,mpileupexample,callsnpsandshortindelsforonediploidindividual:samtoolsmpileup-ufref.faaln.bam|bcftoolsview-bvcg-var.raw.bcfusuallyadd-c50tompileupifmappingqualityisoverestimated,eg.usingbwa-shortbcftoolsviewvar.raw.bcf|vcfutils.plvarfilter-d100var.flt.vcf-dcontrolsthemaxreaddepth,whichshouldbeadjustedtoabouttwicetheaveragereaddepth,mpileupexample,callsnpsandshortindelsformultiplediploidindividuals:samtoolsmpileup-pillumina-ugfref.fa*.bam|bcftoolsview-bvcg-var.raw.bcfindividualsareidentifiedfromthesmtagsinthergheaderlines.-pspecifiesthatindelcandidatesshouldbecollectedonlyfromreadgroupswithrg-pltagsettoilluminabcftoolsviewvar.raw.bcf|vcfutils.plvarfilter-d2000var.flt.vcf,mpileupexample,derivetheafs(allelefrequencyspectrum)onalistofsitesfrommultipleindividuals:samtoolsmpileup-igfref.fa*.bamall.bcfbcftoolsview-blsites.listall.bcfsites.bcfbcftoolsview-cgpcond2sites.bcf/dev/null2sites.1.afsbcftoolsview-cgpsites.1.afssites.bcf/dev/null2sites.2.afsbcftoolsview-cgpsites.2.afssites.bcf/dev/null2sites.3.afs,mpileupexample,generatetheconsensussequence:samtoolsmpileup-ufref.faaln.bam|bcftoolsview-cg-|vcfutils.plvcf2fqcns.fqdumpbaqappliedalignmentforothersnpcallers:samtoolscalmd-braln.bamaln.baq.bam,java版本picard,类似samtools功能的java版本samtools可以在流上（管道）操作，picard好像一般对文件操作。picard的markduplicates能处理samtools的rmdup不能处理的一些问题。（注：markduplicates依赖mate的坐标）,reference,samformat,reference,samtoolsprotocal,bedtools,bed文件格式,bedformat是ucscgenomebrowser和ensemblgenomebrowser都支持的，描述annotationtrack信息的文件格式。它是tab分隔的文本文件，可以有3-12列。前3列是必需的：1chromnameofchr.orscaffold2chromstartstartpos.ofthefeature(include)(0-based)3chromendendpos.ofthefeature(notinclude),bed文件格式,第4-9列是可选的：4namethenameofthebedline5score01000,higher,darker6strand7thickstartstartingpos.drawnthickly8thickendendingpos.drawnthickly9itemrgbdisplaycolor,bed文件格式,第10-12列也是可选的，ucsc使用它们，而ensembl忽略它们：10blockcount#blocks(exons)intheline11blocksizescomma-separated,blocksizes,shouldcorrespondtoblockcount12blockstartscomma-separated,blockstarts,shouldcorrespondtoblockcounteg.snps.bedchr1100101a/g100chr1200201c/g1000,bedgraph文件格式,4列tab分隔的文本文件1chromname2chromstart3chromend4datavalue,bedtools,文件格式转换(bam-bed)bamtobed计算交(overlap)、差、补intersectbedpairtobedpairtopairwindowbed,subtractbedcomplementbed,bedtools,改动bed文件mergebedslopbedshufflebedsortbed计算coveragecoveragebedgenomecoveragebed,bedtools,其他closestbedfastafrombedmaskfastafrombedlinksbed分组简单统计groupby-hdisplayhelppage,bamtobed,bamtobed-iin.bamout.bedconvertbamtobed6format-bed12converttobed12format-bedpeconverttobedpeformat-eduseeditdistance(nmtag)asbedscore,intersectbed,intersectbed-aa.bed-bb.bedreportoverlapsbetweentwobedfiles-cforeacha,report#overlapswithb(count)-ureportoriginalaonceifanyoverlapsfoundinb-vonlyreportthoseathathavenooverlapswithb(similartogrep-v),intersectbed,intersectbed-aa.bed-bb.bed-wareportentire,originalentryinaforeachoverlap-wbreportentire,originalentryinbforeachoverlap,intersectbed,intersectbed-aa.bed-bb.bed-wa-wbreportentire,originalentryinaandbforeachoverlap-woreportoriginalaandband#bpofoverlapbetweenthetwofeatures,intersectbed,intersectbed-aa.bed-bb.bed-ffloatminimumoverlaprequiredasafractionofa-ffloat-rrequirethatthefractionoverlapbereciprocalforaandb-sforcestrandedness,onlyreportwhenaandbonthesamestrand,intersectbed,intersectbed-abama.bam-bb.bedout.baminputfileaisinbam,andouputwillbebam-ubamwriteuncompressedbamoutput-bedwhenusingbaminput,writeoutputasbed-a可跟stdin，表示从stdin读入,intersectbed,eg.findgenesthatoverlaplinesbutnotsines:intersectbed-agenes.bed-blines.bed|intersectbed-astdin-bsines.bed-vretainsebamthatoverlapexons:intersectbed-abamreads.bam-bexons.bedreads.touchingexons.bam,intersectbed,screennovelsnpintersectbed-asnp.calls.bed-bdbsnp.bed-v|intersectbed-astdin-b1kg.bed-vsnp.calls.novel.bed,pairtobed,pairtobed-aa.bedpe-bb.bedreportoverlapsbetweenabedpeandabed-typeeither(default)reportoverlapsifeitherendofaoverlapsb-typeneitherreportaifneitherendofaoverlapsb-sforcestrandedness,pairtobed,pairtobed-aa.bedpe-bb.bed-typebothreportoverlapsifbothendsofaoverlapb-typenotbothreportoverlapsifneitheroronlyoneendofaoverlapb-typexorreportoverlapsifonlyoneendofaoverlapb,pairtobed,pairtobed-aa.bedpe-bb.bed-typeispanreportoverlapsbetweenend1,start2ofa,b-typeospanreportoverlapsbetweenstart1,end2ofa,b-typenotispanreportaifispanofadoesnotoverlapb-typenotospanreportaifospanofadoesnotoverlapb,pairtobed,pairtobed-abama.bam-bb.bedout.baminputfileaispebam,andoutputwillbeinbam-ubamwriteuncompressedbamoutput-bedpewhenusingbaminput,writeoutputasbedpe,pairtobed,eg.returnallstructuralvariantsthatoverlapwithgenesoneitherend:pairtobed-asv.bedpe-bgenes.bedsv.genesretainonlypealignmentswhereneitheroronlyoneendoverlapssegmentalduplications:pairtobed-abamreads.bam-bsegdups.bed-typenotbothreads.notbothssrs.bam,pairtopair,pairtopair-aa.bedpe-bb.bedpereportoverlapsbetweentwobedpe-typeboth(default)reportoverlapsifbothendsofaoverlapb-typeeitherreportoverlapsifeitherendsofaoverlapb-typeneitherreportoverlapsifneitherendofaoverlapsb,pairtopair,eg.findallsvsinsample1thatarealsoinsample2:pairtopair-a1.sv.bedpe-b2.sv.bedpe|cut-f1-101.sv.in2.bedpefindallsvsinsample1thatarenotinsample2:pairtopair-a1.sv.bedpe-b2.sv.bedpe-typeneither|cut-f1-101.sv.notin2.bedpe,windowbed,windowbed-aa.bed-bb.bedexaminesawindowaroundeacha,andreportsallfeaturesinbthatoverlapthewindowforeachoverlap,aandbarereported-wintbpaddedupstreamanddownstreamofeacha(symtericalwindowsarounda),windowbed,windowbed-aa.bed-bb.bed-lintbpaddedupstream(left)ofeacha(assymtericalwindowsarounda)-rintbpaddeddownstream(right)ofeacha(assymtericalwindowsarounda)-swdefine-land-rbasedonstrand,windowbed,windowbed-aa.bed-bb.bed-smonlyreporthitsinbthatoverlapaonthesamestrand-abamfile-ubam-bed,windowbed,eg.reportallsnpsthatarewithin5000bpupstreamor1000bpdownstreamofgenes:windowbed-agenes.bed-bsnps.bed-l5000-r1000-sw,subtractbed,subtractbed-aa.bed-bb.bedremovesathatoverlapbyanyfeatureinb-ffloatminimumoverlaprequiredasafractionofa-sforcestrandednesseg.removeintronsfromgenefeatures:subtractbed-agenes.bed-bintron.bed,complementbed,complementbed-iin.bed-ggenomereturnbpcomplementofthefeaturefilegenomefileshouldtab-delimitedandcontain:chromnamechromsizehowtogetchromsize?see-heg.reportallintervalsinhumangenomethatarenotcoveredbyrepetitiveelements:complementbed-irepeatmasker.bed-ghg18.genome,mergebed,mergebed-iin.bedmergeoverlappingentriesintoasingleentry-nreport#bedentriesmerged1isreportedifnomergingoccurred-dintmaximumdistance(bp)betweenfeaturesallowedtobemerged0,mergebed,mergebed-iin.bed-sforcestrandedness-scoresstrhowtodeterminethescoresofthemergedsum,min,max,mean,median,mode,antimode,collapse,slopbed,slopbed-iin.bed-ggenome将feature向两侧扩展-bint/floatineachdirection-lint/floatinleft(upstream)direction-rint/floatinright(downstream)direction,slopbed,slopbed-iin.bed-ggenome-sdefine-land-rbasedonstrand-pctspecifyafractionofthefeatureslength,insteadofabsolutebasesonlywith-pct,-b/-l/-rcanbegivenfloateg.slopbed-iprobes.bed-ghg19.genome-b500,shufflebed,shufflebed-iin.bed-ggenomerandomlypermutethelocationsoffeaturesamongagenome-exclbed_filepreventplacingintheseintervals-ffloatmaximumoverlapfractionofthefeaturewithan-exclfeature-chromkeepfeaturesonthesamechromosome,sortbed,sortbed-iin.bedsortbedfile,bychrom,thenstartposition输出排序好的bed到stdoutsee-hformoreoptions,coveragebed,coveragebed-aa.bed-bb.bedcomputethedepthandbreadthofcoverageoffeaturesfromaonintervalsinbdefaultoutput:(foreachentryinb)#aoverlappedthebinterval#basesinbhad0coveragethelengthofbthefractionofbasesinbhad0coverage,coveragebed,coveragebed-aa.bed-bb.bed-sforcestrandedness-dreportthedepthateachposition(base)ineachbputecoverageandcreateabedgraphcoveragebed-areads.bed-bwin10kb.bed|cut-f1-4win10kb.cov.bedg,coveragebed,bydefault,coveragebedcountsanyfeatureinathatoverlapsbby=1bp.ifyouwanttospecifyaminimumfraction,youcanfirstuseintersectbed:intersectbed-aa.bed-bb.bed-f1.0|coveragebed-astdin-bb.beda_totally_in_b.coverage,coveragebed,workwithsamtoolsorbamfilebamtobed-ireads.bam|coveragebed-astdin-bexons.bedexons.bed.coveragesamtoolsview-bf0 x2reads.bam|bamtobed-istdin|coveragebed-astdin-bexons.bedexons.bed.pair_proper_mapped.coverage,coveragebed,computeseparatedlyforeachstrandbamtobed-ireads.bam|grep+$|coveragebed-astdin-bgenes.bedgenes.bed.forward.coveragebamtobed-ireads.bam|grep-$|coveragebed-astdin-bgenes.bedgenes.bed.reverse.coverage,genomecoveragebed,genomecoveragebed-iin.bed-ggenomecomputethecoverageamongagenomeinputbedfilemustbegroupedbychr.defaultoutput:(foreachchromosomeinb)depth(0max)#basesinthechr.hadthisdepththelengthofthechr.thefractionofbasesinthechr.thathadthisdepth-strand+/-,genomecoveragebed,-ibamfileinputinbammustbesortedbyposition(samtoolssort)-dreportthedepthateachposition(base)-bgreportdepthinbedgraphformat-bgareportdepthinbedgraphformat,alsoreportzerocoverageregions,genomecoveragebed,eg.genomecoveragebed-ibambwa_result.sorted.bam-ghuman_chrm.genomegenomecoveragebed-ibambwa_result.sorted.bam-ghuman_chrm.genome|grepchrm|awkbeginsum=0sum+=($2*$5)endprintsum,closestbed,closestbed-aa.bed-bb.bedforeacha,findstheclosestfeatureinb-tall(default)reportallties-tfirstreportthefirsttiethatoccurredinb-tlastreportthelasttiethatoccurredinb,closestbed,closestbed-aa.bed-bb.bed-sforcestrandedness,findtheclosestfeatureonthesamestrandeg.findtheclosestalutoeachgene:closestbed-agenes.bed-balus.bed,fastafrombed,maskfastafrombed-fiin.fasta-bedextract.bed-foout.fasta根据extract.bed文件，从in.fasta中抽取一些区域，输出到out.fasta-sforcestrandedness-nameusenamefieldforfastaheader,maskfastafrombed,maskfastafrombed-fiin.fasta-bedmask.bed-foout.fasta根据mask.bed文件，mask掉in.fasta中的一些区域，输出到out.fasta-