【6.1.2.1】samtools
对于一个生物信息从业人员来说,玩的就是序列,这里面离不开的是序列的比对,作为一个强大的查看序列比对结果的工具,有必要好好认识一下samtools。
一、samtools的简介
SAMtools: Primer / Tutorial http://biobits.org/samtools_primer.html
1.1 samtools的安装
下载samtools工具:http://sourceforge.net/projects/samtools/files/
wget -c https://nchc.dl.sourceforge.net/project/samtools/samtools/1.9/samtools-1.9.tar.bz2
tar -xjf samtools-1.9.tar.bz2
cd samtools-1.9/
./configure
make
make install
#samtoolpath=`pwd`
#PATH=PATH:$samtoolpath
修改bash文件:
export SAMTOOLS_HOME=/Users/ecerami/libraries/samtools-1.9
export PATH=$SAMTOOLS_HOME:$PATH
export PATH=$SAMTOOLS_HOME/bcftools/:$PATH
export PATH=$SAMTOOLS_HOME/misc/:$PATH
二、常用命令
2.1 help
Command: view SAM<->BAM conversion 两种文件的互换
sort sort alignment file
mpileup multi-way pileup
depth compute the depth
faidx index/extract FASTA
tview text alignment viewer
index index alignment
idxstats BAM index stats (r595 or later)
fixmate fix mate information
flagstat simple stats
calmd recalculate MD/NM tags and '=' bases
merge merge sorted alignments
rmdup remove PCR duplicates
reheader replace BAM header
cat 只支持合并BAM文件
bedcov read depth per BED region
targetcut cut fosmid regions (for fosmid pool only)
phase phase heterozygotes
bamshuf shuffle and group alignments by name
2.2 View
将sam文件转换成bam文件;然后对bam文件进行各种操作,比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的,因而当输入为sam文件的时候,不能进行该操作);最后将排序或提取得到的数据输出为bam或sam(默认的)格式。
bam文件优点:bam文件为二进制文件,占用的磁盘空间比sam文本文件小;利用bam二进制文件的运算速度快。
Usage: samtools view [options] | [region1 [...]]
默认情况下不加 region,则是输出所有的 region.
Options: -b output BAM
默认下输出是 SAM 格式文件,该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header,该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
-S input is SAM 。 默认下输入是 BAM 文件,若是输入是 SAM 文件,则最好加该参数,否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数,能节约时间,但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ ,
are taken into account.
-1 fast compression (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-s FLOAT fraction of templates to subsample; integer part as seed [-1]
-? longer help
samtools 将SAM文件转化为BAM文件
samtools view -bS eg2.sam > eg2.bam
也可以输入:
samtools view -b -S -o alignments/sim_reads_aligned.bam alignments/sim_reads_aligned.sam
-b: indicates that the output is BAM.
-S: indicates that the input is SAM.
-o: specifies the name of the output file.
因为bam文件是一个压缩的二进制文件,不能直接看,所以 需要使用view功能
samtools view alignments/sim_reads_aligned.bam | more
也可以使用view来仅仅展示那些匹配到specific filtering criteria的reads
samtools view -f 4 alignments/sim_reads_aligned.bam | more
-f int :仅仅提取匹配specified SAM flag的reads,上面的那个例子中,我们仅仅提取那些flag value 为4的匹配上去的reads
samtools view -F 4 alignments/sim_reads_aligned.bam | more
-F init: 上面的例子中是删掉匹配上flag value 为4的reads
samtools view -c -f 4 alignments/sim_reads_aligned.bam
上面不产生结果文件,仅仅甘肃你多少个reads没有匹配到reference genome
samtools view -q 42 -c alignments/sim_reads_aligned.bam
显示匹配的分数在42分以上的reads数
#提取比对到参考序列上的比对结果
$ samtools view -bF 4 abc.bam > abc.F.bam
提取paired
reads中两条reads都比对到参考序列上的比对结果,只需要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 abc.bam > abc.F12.bam
#提取没有比对到参考序列上的比对结果
samtools view -bf 4 abc.bam > abc.f.bam
#提取bam文件中比对到caffold1上的比对结果,并保存到sam文件格式
samtools view abc.bam scaffold1 > scaffold1.sam
#提取scaffold1上能比对到30k到100k区域的比对结果
samtools view abc.bam scaffold1:30000-100000 > scaffold1_30k-100k.sam
#根据fasta文件,将 header 加入到 sam 或 bam 文件中
samtools view -T genome.fasta -h scaffold1.sam > scaffold1.h.sam
This is not a direct answer to your question, but you can do some basic counts using samtools
注意f大小写的区分,以及后面跟的数字0
number of mapped reads:
(既可以读sam文件,也可以读bam文件)
samtools view -S -F0x4 -c reads.sam >mapped_reads
number of unmapped reads:
samtools view -S -f0x4 -c reads.sam
number of uniquely mapped reads:
samtools view -S -F0x4 -c -q 1 reads.sam
map到某个bed文件的reads
samtools view -b -L brca_LP.target.zerobased.bed S1570169_sort.bam >aa.bam
samtools view -S -F0x4 -c aa.bam>mapped_reads
2.3 sort,index
用samtools排序和然后建立索引
samtools sort eg2.bam eg2.sorted
samtools index eg2.sorted.bam
第一步生成eg2.sorted.bam文件,第二步生成sim_reads_aligned.sorted.bam.bai,如果没有排序的话,是没有办法建立索引的。 建立索引后将产生后缀为.bai的文件,用于快速的随机处理。很多情况下需要有bai文件的存在,特别是显示序列比对情况下。比如samtool的tview命令就需要;gbrowse2显示reads的比对图形的时候也需要。
2.4 merge/cat
merge将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件不需要则是已经sort过了的。
Usage: samtools merge [-nr] [-h inh.sam] [...]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to [in1.bam]
Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users
must provide the correct header with -h, or uses Picard which properly maintains
the header dictionary in merging.
cat 连接多个bam文件,适用于非sorted的bam文件
$ samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
2.5 faidx
对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条(子)序列
Usage: samtools faidx [ [...]]
对基因组文件建立索引
$ samtools faidx genome.fasta
#生成了索引文件genome.fasta.fai,是一个文本文件,分成了5列。第一列是子序列的名称;第二列是子序列的长度;个人认为“第三列是序列所在的位置”,因为该数字从上往下逐渐变大,最后的数字是genome.fasta文件的大小;第4和5列不知是啥意思。于是通过此文件,可以定位子序列在fasta文件在磁盘上的存放位置,直接快速调出子序列。
#由于有索引文件,可以使用以下命令很快从基因组中提取到fasta格式的子序列
$ samtools faidx genome.fasta scffold_10 > scaffold_10.fasta
2.6 tview
tview能直观的显示出reads比对基因组的情况,和基因组浏览器有点类似。
Usage: samtools tview [ref.fasta]
当给出参考基因组的时候,会在第一排显示参考基因组的序列,否则,第一排全用N表示。
按下 g ,则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第
10号scaffold的第1000个碱基位点处。
使用H(左)J(上)K(下)L(右)移动显示界面。大写字母移动快,小写字母移动慢。
使用空格建向左快速移动(和 L 类似),使用Backspace键向左快速移动(和 H 类似)。
Ctrl+H 向左移动1kb碱基距离; Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量,碱基质量,核苷酸等。30~40的碱基质量或比对质量使用白色表示;
20~30黄色;10~20绿色;0~10蓝色。
使用点号'.'切换显示碱基和点号;使用r切换显示read name等
还有很多其它的使用说明,具体按 ? 键来查看。
2.7 flagstat 统计比对结果
samtools flagstat test.sam > test.flagstat
说明:因为使用的bwa的mem比对方法生成的bam文件,里面会保留sencondly比对的结果。所以使用flagstat给出的结果会有偏差。
flagstat使用方法:samtools flagstat *.bam>*.stat
输出结果:
20607872 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数
0 + 0 duplicates
19372694 + 0 mapped (94.01%:-nan%) ##总体上reads的匹配率
20607872 + 0 paired in sequencing #有多少reads是属于paired reads
10301155 + 0 read1 ##reads1中的reads数
10306717 + 0 read2 #falg为2出现的次数;
11228982 + 0 properly paired (54.49%:-nan%)
#完美匹配的reads数:比对到同一条参考序列,并且两条reads之间的距离符合设置的阈值
18965125 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
407569 + 0 singletons (1.98%:-nan%)
#单独一条匹配到参考序列上的reads数,和上一个相加,则是总的匹配上的reads数。
3059705 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
1712129 + 0 with mate mapped to a different chr (mapQ>=5)
2.8 depth 得到每个碱基位点的测序深度
得到每个碱基位点的测序深度,并输出到标准输出。
samtools index test_sort.bam test_sort.bai
samtools depth -b /home/qqin/genopipe_py/genopipe/data/BRCAWise.exons.bed test_sort.bam > test_exons.depth
awk 'BEGIN{OFS="\t"}{if($3>20) print $1,$2,$2,$3}' test.depth > test.cut20.depth
得到单个碱基覆盖度为20的碱基
bedtools merge -i test.cut20.depth > test.cut20.merged.region
awk 'BEGIN{OFS="\t"}{print $1,$2+1,$3-1 > "test.cut20.merged.region"}'
test.cut20.merged.region
samtools depth -b /home/qqin/genopipe_py/genopipe/data/BRCAWise.exons.bed test_sort.bam -m 11100000000 > test_exons.depth
如果测序深度太深,需要制定-m
2.9. reheader 替换bam文件的头
$ samtools reheader <in.header.sam> <in.bam>
2.10 idxstats 评估map的质量
samtools idxstats aln.sorted.bam
注意,这一步之前需要经过sort和index。结果会显示:
chr1 195471971 6112404 0
chr10 130694993 3933316 0
chr11 122082543 6550325 0
chr12 120129022 3876527 0
chr13 120421639 5511799 0
chr14 124902244 3949332 0
chr15 104043685 3872649 0
chr16 98207768 6038669 0
chr17 94987271 13544866 0
chr18 90702639 4739331 0
chr19 61431566 2706779 0
chr2 182113224 8517357 0
chr3 160039680 5647950 0
chr4 156508116 4880584 0
chr5 151834684 6134814 0
chr6 149736546 7955095 0
chr7 145441459 5463859 0
chr8 129401213 5216734 0
chr9 124595110 7122219 0
chrM 16299 1091260 0
chrX 171031299 3248378 0
chrY 91744698 259078 0
* 0 0 0
其中第一列是染色体名称,第二列是序列长度,第三列是mapped reads数,第四列是unmapped reads数,paired reads中有一端能匹配到该scaffold上,而另外一端不匹配到任何scaffolds上的reads数。idxstats 统计一个表格,4列,分别为”序列名,序列长度,比对上的reads数,unmapped reads number”。第4列应该是paired reads中有一端能匹配到该scaffold上,而另外一端不匹配到任何scaffolds上的reads数。
2.11 mpileup 寻找突变通过VCF格式
samtools还有个非常重要的命令mpileup,以前为pileup。该命令用于生成bcf文件,再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件,在samtools的安装文件夹中可以找到。
最常用的参数有2:
-f 来输入有索引文件的fasta参考序列;
-g 输出到bcf格式。用法和最简单的例子如下
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
$ samtools mpileup -f genome.fasta abc.bam > abc.txt
$ samtools mpileup -gSDf genome.fasta abc.bam > abc.bcf
$ samtools mpileup -guSDf genome.fasta abc.bam | \
bcftools view -cvNg - > abc.vcf
mpileup不使用-u或-g参数时,则不生成二进制的bcf文件,而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况;该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如:
scaffold_1 2841 A 11 ,,,...,.... BHIGDGIJ?FF
scaffold_1 2842 C 12 ,$,,...,....^I. CFGEGEGGCFF+
scaffold_1 2843 G 11 ,,...,..... FDDDDCD?DD+
scaffold_1 2844 G 11 ,,...,..... FA?AAAA<AA+
scaffold_1 2845 G 11 ,,...,..... F656666166*
scaffold_1 2846 A 11 ,,...,..... (1.1111)11*
scaffold_1 2847 A 11 ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG %.+....-..)
scaffold_1 2848 N 11 agGGGgGGGGG !!$!!!!!!!!
scaffold_1 2849 A 11 c$,...,..... !0000000000
scaffold_1 2850 A 10 ,...,..... 353333333
mpileup生成的结果包含6行:参考序列名;位置;参考碱基;比对上的reads数;比对情况;比对上的碱基的质量。其中第5列比较复杂,解释如下:
-
‘.’代表与参考序列正链匹配。
-
‘,’代表与参考序列负链匹配。
-
‘ATCGN’代表在正链上的不匹配。
-
‘atcgn’代表在负链上的不匹配。
-
‘*’代表模糊碱基
-
‘^’代表匹配的碱基是一个read的开始;’^‘后面紧跟的ascii码减去33代表比对质量;这两个符号修饰的是后面的碱基,其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。
-
‘$’代表一个read的结束,该符号修饰的是其前面的碱基。
-
正则式’+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基;比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。
-
正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基;
samtools mpileup -uf $BT2_HOME/example/reference/lambda_virus.fa eg2.sorted.bam | bcftools view -bvcg - > eg2.raw.bcf
samtools mpileup -g -f genomes/NC_008253.fna alignments/sim_reads_aligned.sorted.bam > variants/sim_variants.bcf
通过工具mpileup统计通过比对基因型相似的reads
-g: 指导SAMtools产生基因型相似性的二进制访问格式(binary call format)(BCF).
-f: 指导SAMtools使用指定的相关基因组,相关基因组必须被指定,上面的例子指定的是E. coli
mpileup命令能够自动浏览咩个位置的比对的reads,计算所有基因型的可能性,然后统计在我们样品中这些基因型的每一个的可能性。
2.12 bcftools突变位点
bcftools和samtools类似,用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件,后者为其二进制文件。 bcftools使用简单,最主要的命令是view命令,其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为:
$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]
[-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]
[-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]
[-T trioType] in.bcf [region]
$ bcftools view -cvNg abc.bcf > snp_indel.vcf
生成的结果文件为vcf格式,有10列,分别是:
1 参考序列名;
2 variant所在的left-most位置;
3 variant的ID(默认未设置,用’.'表示);
4 参考序列的allele;
5 variant的allele(有多个alleles,则用’,'分隔);
6 variant/reference QUALity;
7 FILTers applied;
8 variant的信息,使用分号隔开;
9 FORMAT of the genotype fields, separated by colon (optional);
10 SAMPLE genotypes and per-sample information (optional)。
例如:
scaffold_1 2847 . A AACGGTGAAG 194 . INDEL;DP=11;VDB=0.0401;AF1=1;AC1=2;DP4=0,0,8,3;MQ=35;FQ=-67.5 GT:PL:GQ 1/1:235,33,0:63
scaffold_1 3908 . G A 111 . DP=13;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,5,7;MQ=42;FQ=-63 GT:PL:GQ 1/1:144,36,0:69
scaffold_1 4500 . A G 31.5 . DP=8;VDB=0.0034;AF1=1;AC1=2;DP4=0,0,1,3;MQ=42;FQ=-39 GT:PL:GQ 1/1:64,12,0:21
scaffold_1 4581 . TGGNGG TGG 145 . INDEL;DP=8;VDB=0.0308;AF1=1;AC1=2;DP4=0,0,0,8;MQ=42;FQ=-58.5 GT:PL:GQ 1/1:186,24,0:45
scaffold_1 4644 . G A 195 . DP=21;VDB=0.0198;AF1=1;AC1=2;DP4=0,0,10,10;MQ=42;FQ=-87 GT:PL:GQ 1/1:228,60,0:99
scaffold_1 4827 . NACAAAGA NA 4.42 . INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=40;FQ=-37.5 GT:PL:GQ 0/1:40,3,0:3
scaffold_1 4854 . A G 48 . DP=6;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,2,1;MQ=41;FQ=-36 GT:PL:GQ 1/1:80,9,0:16
scaffold_1 5120 . A G 85 . DP=8;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,5,3;MQ=42;FQ=-51 GT:PL:GQ 1/1:118,24,0:45
第8列中显示了对variants的信息描述,比较重要,其中的 Tag 的描述如下:
Tag Format Description
AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
DP int Raw read depth (without quality filtering)
DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
MQ int Root-Mean-Square mapping quality of covering reads
PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
QCHI2 int Phred-scaled PCHI2
RP int # permutations yielding a smaller PCHI2
CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
UGT string Most probable genotype configuration without the trio constraint
CGT string Most probable configuration with the trio constraint
bcftools view 的具体参数如下:
Input/Output Options:
-A Retain all possible alternate alleles at variant sites. By default, the view command discards unlikely alleles.
-b Output in the BCF format. The default is VCF.
-D FILE Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
-F Indicate PL is generated by r921 or before (ordering is different).
-G Suppress all individual genotype information.
-l FILE List of sites at which information are outputted [all sites]
-N Skip sites where the REF field is not A/C/G/T
-Q Output the QCALL likelihood format
-s FILE List of samples to use. The first column in the input gives the sample names and the second gives the ploidy, which can only be 1 or 2. When the 2nd column is absent, the sample ploidy is assumed to be 2. In the output, the ordering of samples will be identical to the one in FILE. [null]
-S The input is VCF instead of BCF.
-u Uncompressed BCF output (force -b).
Consensus/Variant Calling Options:
-c Call variants using Bayesian inference. This option automatically invokes option -e.
-d FLOAT When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
当有多个sample用于variants calling时,比如多个转录组数据或多个重测序
数据需要比对到参考基因组上,设置该值,表明至少有该比例的
samples在该位点都有覆盖才计算入variant.所以对于只有一个sample的情况
下,该值设置在0~1之间没有意义,大于1则得不到任何结果。
-e Perform max-likelihood inference only, including estimating the site allele frequency, testing Hardy-Weinberg equlibrium and testing associations with LRT.
-g Call per-sample genotypes at variant sites (force -c)
-i FLOAT Ratio of INDEL-to-SNP mutation rate [0.15]
-p FLOAT A site is considered to be a variant if P(ref|D)
-t FLOAT Scaled muttion rate for variant calling [0.001]
-T STR Enable pair/trio calling. For trio calling, option -s is usually needed to be applied to configure the trio members and their ordering. In the file supplied to the option -s, the first sample must be the child, the second the father and the third the mother. The valid values of STR are ‘pair’, ‘trioauto’, ‘trioxd’ and ‘trioxs’, where ‘pair’ calls differences between two input samples, and ‘trioxd’ (‘trioxs’) specifies that the input is from the X chromosome non-PAR regions and the child is a female (male). [null]
-v Output variant sites only (force -c)
Contrast Calling and Association Test Options:
-1 INT Number of group-1 samples. This option is used for dividing the samples into two groups for contrast SNP calling or association test. When this option is in use, the following VCF INFO will be outputted: PC2, PCHI2 and QCHI2. [0]
-U INT Number of permutations for association test (effective only with -1) [0]
-X FLOAT Only perform permutations for P(chi^2)
使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行尤为重要,提供了4个数据:1 比对结果和正链一致的reads数、2 比对结果和负链一致的reads数、3 比对结果在正链的variant上的reads数、4 比对结果在负链的variant上的reads数。可以设定 (value3 + value4)大于某一阈值,才算是variant。比如:
$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ && ($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf > snp_indel.final.vcf
bcftools view eg2.raw.bcf
bcftools view -c -v variants/sim_variants.bcf > variants/sim_variants.vcf
-c: directs bcftools to call SNPs.
-v: directs bcftools to only output potential variants
bcftools view命令用之前产生的基因型相似性来寻找SNP和节点,结果保存在vcf中
2.13 rmdup
NGS上机测序前需要进行PCR一步,使一个模板扩增出一簇,从而在上机测序的时候表现出为1个点,即一个reads。若一个模板扩增出了多簇,结果得到了多个reads,这些reads的坐标(coordinates)是相近的。在进行了reads比对后需要将这些由PCR duplicates获得的reads去掉,并只保留最高比对质量的read。使用rmdup命令即可完成.
Usage: samtools rmdup [-sS]
-s 对single-end reads。默认情况下,只对paired-end reads
-S 将Paired-end reads作为single-end reads处理。
$ samtools input.sorted.bam output.bam
三、讨论
3.1 depth和mpileup的区别
depth: " compute the per-base depth"
mpileup: "the number of reads covering the site"
被双链覆盖
被单链覆盖
mpileup
-A, --count-orphans do not discard anomalous read pairs
discard anomalous read pairs (samtools mpileup -A)
3.2 samtools和igv看到的depth不一致
这里面有两个原因:序列的质量以及duplicates
samtools在计算depth的时候,会过滤掉duplicates的序列,解决办法:
bedtools genomecov -ibam GWB0514095_sort.primerclipped.realign.recalibrate.bam -bg |grep 4119781
17 41197808 41197811 5802
17 41197811 41197812 5800
17 41197812 41197821 5802
因为得到的是真实的bed文件,实际上是17:41197809-41197811:5802
3.3 将bam文件转换为fastq文件
有时候,我们需要提取出比对到一段参考序列的reads,进行小范围的分析,以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序:http://www.hudsonalpha.org/gsl/information/software/bam2fastq
$ wget http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz
$ tar zxf bam2fastq-1.1.0.tgz
$ cd bam2fastq-1.1.0
$ make
$ ./bam2fastq <in.bam>
3.4 报错: hts_open_format
samtools sort NA18562.mapped.ILLUMINA.bwa.CHB.low_coverage.20101123.bam -o NA18562.mapped.
ILLUMINA.bwa.CHB.low_coverage.20101123.sorted.bam -m 4G --threads 20
[E::hts_open_format] fail to open file 'NA18562.mapped.ILLUMINA.bwa.CHB.low_coverage.20101123.bam.tmp.0000.bam'
问题原因:
I have had this exact same problem recently. It turned out to be that the system limits the number of opened files. After increasing the limit, the error is gone.
解决办法:
- 去1000g数据库中下载对应的bai文件(。。。。)
- 修改系统容许打开的文件,提高命令中threads的值(未这么做)
四、我的案例
3.1 求得coverage
cd /sam/uncltured/contig7/index #建立索引
/sam/bowtie2/bowtie2-build ../contig.fasta scaffold #这里的contig.fasta是我变形以后的scaffold,
先用paired 去比对
/sam/bowtie2/bowtie2 -p 12 -x scaffold -1 /sam/idba_ud/myproject/reads2/output_forward_paired.fq -2 /sam/idba_ud/myproject/reads2/output_reverse_paired.fq -S eg2.sam
219079809 reads; of these:
219079809 (100.00%) were paired; of these:
9992084 (4.56%) aligned concordantly 0 times
190538854 (86.97%) aligned concordantly exactly 1 time
18548871 (8.47%) aligned concordantly >1 times
---
9992084 pairs aligned concordantly 0 times; of these:
985097 (9.86%) aligned discordantly 1 time
----
9006987 pairs aligned 0 times concordantly or discordantly; of these:
18013974 mates make up the pairs; of these:
9341217 (51.86%) aligned 0 times
4121203 (22.88%) aligned exactly 1 time
4551554 (25.27%) aligned >1 times
97.87% overall alignment rate
再用unpair去比对
/sam/bowtie2/bowtie2 -p 12 -x scaffold -U /sam/idba_ud/myproject/reads2/unpair.fq -S eg1.sam
9910673 reads; of these:
9910673 (100.00%) were unpaired; of these:
303363 (3.06%) aligned 0 times
8085141 (81.58%) aligned exactly 1 time
1522169 (15.36%) aligned >1 times
96.94% overall alignment rate
samtools view -bS eg1.sam > eg1.bam
[samopen] SAM header is present: 144759 sequences.
samtools view -bS eg2.sam > eg2.bam [11:23PM]
[samopen] SAM header is present: 144759 sequences
samtools cat eg1.bam eg2.bam >eg.bam
samtools sort eg.bam eg.sorted
[bam_sort_core] merging from 184 files...
参考资料
- (推荐)SAMtools: Primer / Tutorial http://biobits.org/samtools_primer.html
- 菜菜的一天的博客 http://blog.sina.com.cn/s/blog_670445240101l30k.html
- 糗世界的博客
- https://www.biostars.org/p/57968/
- https://www.plob.org/article/7112.html
- https://mp.weixin.qq.com/s?__biz=MzI2MjA1MDQxMg==&mid=2649709965&idx=2&sn=a5a7b2de70ad0fac61eb541946a1a8f9&chksm=f24afa8ec53d739824702496dc642693c05c2af40a85a590921c46ee403c14d90b0abdd0ef67&scene=21#wechat_redirect
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn