必須學(xué)習(xí)這篇英文：http://www.htslib.org/doc/samtools.html
【繼續(xù)更新ing】
samtools view -h
查看bam文件囊嘉，包含頭文件蹭秋，去掉-h菱父，不包含
samtools view -h s.bam|less -S
samtools view s.bam|less -S
# 提取chr1染色體脯丝，生成只有chr1的bam文件
samtools view -h -b s.bam chr1 >s.chr1.bam
samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
samtools view -c
-c print only the count of matching records；統(tǒng)計(jì)比對(duì)條數(shù);結(jié)果等于samtools flagstat的總reads條數(shù)雳灵。
# 單端比對(duì)reads數(shù)目的統(tǒng)計(jì)
cuiqingmei 2019/10/09 10:57:10 /ifs9/BC_PS/cuiqingmei/ass/3.mapping/result_bam
$ samtools view -c CL100141207_L01_69.sort.bam
31504107
cuiqingmei 2019/10/09 10:57:38 /ifs9/BC_PS/cuiqingmei/ass/3.mapping/result_bam
$ zcat ../CL100141207_L01_69.fq.gz |wc -l
126016428
cuiqingmei 2019/10/09 11:02:26 /ifs9/BC_PS/cuiqingmei/ass/3.mapping/result_bam
$ echo "scale=4;126016428/4"|bc
31504107.0000
qmcui 12:07:32 /teach/project/1.rna/4.clean_data_25000reads
$ echo `zcat SRR1039510_1_val_1.100000.fq.gz|wc -l`/4|bc
25000
qmcui 12:08:12 /teach/project/1.rna/4.clean_data_25000reads
$ echo `zcat SRR1039510_2_val_2.100000.fq.gz|wc -l`/4|bc
25000
qmcui 12:45:49 /teach/project/1.rna/4.clean_data_25000reads
$ samtools view -c ../5.sort.bam/SRR1039510.sort.bam
55621
“Instead of printing the alignments, only count them and print the total number.” samtools view -c 算的不是reads祈餐，而是alignments數(shù)。當(dāng)read有很多位置可以align上同時(shí)又都輸出了扳还，用samtools view -c 會(huì)比實(shí)際reads樹(shù)木要多～～～ by: 嚴(yán)云
samtools view -F/-f 4
看下表才避，4是unmapped；所以我們-F（none of） 4就是找到map的reads數(shù)目
-q INT only include reads with mapping quality >= INT [0]
-m INT only include reads with number of CIGAR operations consuming
query sequence >= INT [0]
-f INT only include reads with all of the FLAGs in INT present [0]
-F INT only include reads with none of the FLAGS in INT present [0]
-G INT only EXCLUDE reads with all of the FLAGs in INT present [0]
So-f 4 only output alignments that areunmapped(flag 0x0004 is set) and -F 4only output alignments that are not unmapped (i.e. flag 0x0004 is not set), hence these would onlyinclude mapped alignments.
cuiqingmei 2019/10/09 11:02:58 /ifs9/BC_PS/cuiqingmei/ass/3.mapping/result_bam
$ samtools view -c CL100141207_L01_69.sort.bam -F 4
30441976
# 想知道有多少paired end reads有mate并且都有map時(shí)氨距，可以使用-f 1 -F 12來(lái)過(guò)濾桑逝。
samtools view -c -f 1 -F 12 test.bam
# 12是 read unmapped、 mate unmapped;12=8+4
# Mapped reads only
$ samtools view -c -F 4 HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam
5068340
# Unmapped reads only
$ samtools view -c -f 4 HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam
149982
| Flag | Chr | Description |
|---|---|---|
| 0×0001 | p | the read is paired in sequencing |
| 0×0002 | P | the read is mapped in a proper pair |
| 0×0004 | u | the query sequence itself is unmapped |
| 0×0008 | U | the mate is unmapped |
| 0×0010 | r | strand of the query (1 for reverse) |
| 0×0020 | R | strand of the mate |
| 0×0040 | 1 | the read is the first read in a pair |
| 0×0080 | 2 | the read is the second read in a pair |
| 0×0100 | s | the alignment is not primary |
| 0×0200 | f | the read fails platform/vendor quality checks |
| 0×0400 | d | the read is either a PCR or an optical duplicate |
samtools sort
bam文件必須排序
一般按照默認(rèn)參數(shù)俏让，控制線程數(shù)均可楞遏，-@ 5即5個(gè)線程茬暇。
-n Sort by read name：按照reads名字來(lái)排序
-l INT Set compression level, from 0 (uncompressed) to 9 (best)：設(shè)置壓縮倍數(shù)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]：
Usage: samtools sort [-n] [-m <maxMem>] <in.bam> <out.prefix>
-m 參數(shù)默認(rèn)下是 500,000,000 即500M（不支持K，M寡喝，G等縮寫(xiě)）糙俗。對(duì)于處理大數(shù)據(jù)時(shí)，如果內(nèi)存夠用预鬓，則設(shè)置大點(diǎn)的值巧骚，以節(jié)約時(shí)間。
-n 設(shè)定排序方式按short reads的ID排序珊皿。默認(rèn)下是按序列在fasta文件中的順序（即header）和序列從左往右的位點(diǎn)排序。
$ samtools sort -@ 5 -o output.sort.bam input.sam
$ samtools sort -@ 5 -o output.sort.bam input.bam
# 排序sam巨税，bam均可蟋定，而且排序后結(jié)果默認(rèn)生成bam
# samtools sort -T /tmp/aln.sorted -o aln.sorted.bam aln.bam
samtools flags
explain BAM flags：解釋bam文件第二列標(biāo)簽的含義
cuiqingmei 2019/10/09 11:51:33
$ samtools flags 4
0x4 4 UNMAP
cuiqingmei 2019/10/09 11:52:44
$ samtools flags 355
0x163 355 PAIRED,PROPER_PAIR,MREVERSE,READ1,SECONDARY
samtools index
給sam或者bam文件建立索引，一般下游程序需要位置信息的時(shí)候草添，就必須在bam同目錄下存在bai索引文件驶兜。
# samtools 1.8版本
$ samtools index
Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
-b Generate BAI-format index for BAM files [default]
-c Generate CSI-format index for BAM files
-m INT Set minimum interval size for CSI indices to 2^INT [14]
-@ INT Sets the number of threads [none]
cuiqingmei 2019/10/09 11:56:21 /ifs9/
$ samtools index CL100141207_L01_69.sort.bam
samtools idxstat
結(jié)果解釋?zhuān)簉eference sequence name, sequence length, # mapped reads and # unmapped reads，\t分割
$ samtools idxstats CL100141207_L01_69.sort.bam
chr1 249250621 2443169 0
chr2 243199373 2575343 0
chr3 198022430 2018108 0
chr4 191154276 1987277 0
chr5 180915260 1838927 0
...
* 0 0 1062131
samtools markdup
去重
EXAMPLE
# The first sort can be omitted if the file is already name ordered
samtools sort -n -o namesort.bam example.bam
# Add ms and MC tags for markdup to use later
samtools fixmate -m namesort.bam fixmate.bam
# -m:Add ms (mate score) tags. These are used by markdup to select the best reads to keep.
# Markdup needs position order
samtools sort -o positionsort.bam fixmate.bam
# Finally mark duplicates
samtools markdup positionsort.bam markdup.bam
samtools rmdup：過(guò)時(shí)
samtools rmdup [-sS] <input.srt.bam> <out.bam>
This command is obsolete. Use markdup instead.命令過(guò)時(shí)远寸，最好不要再使用抄淑。
值得學(xué)習(xí)翻譯鏈接：
http://qnot.org/2012/04/14/counting-the-number-of-reads-in-a-bam-file/
http://www.htslib.org/doc/samtools.html
還可以學(xué)習(xí)：http://www.chenlianfu.com/?p=1399
