samtools

学习资料

• Samtools

sam转bam

由于sam格式的文本文件过大，一般将其转换为bam格式的文件

samtools view -bS  -@ 10 EV_str.sam -o EV_str.bam

bam文件排序

samtools  sort F11162.bam -@ 10 -o   F11162.sort.bam 

• -@, --threads:要使用的附加线程数

建立bam文件的索引

bam文件必须进行sort后才能进行index

samtools index EV_3.bam

统计比对信息

samtools  flagstat -@ 10  hisat2.bam

20607872 + 0 in total (QC-passed reads + QC-failed reads)   #reads总数。
0 + 0 duplicates     #PCR duplication reads总数。
19372694 + 0 mapped (94.01%:-nan%)    #reads比对率。
20607872 + 0 paired in sequencing   #属于PE read的reads总数。
10301155 + 0 read1             #PE read中Read 1 的reads 总数。
10306717 + 0 read2        #PE read中Read 2 的reads 总数。
11228982 + 0 properly paired (54.49%:-nan%)    完美比对的reads总数。PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值。
18965125 + 0 with itself and mate mapped       #PE两端reads都比对上参考序列的reads总数。
407569 + 0 singletons (1.98%:-nan%)    #PE两端reads，其中一端比上，另一端没比上的reads总数。
3059705 + 0 with mate mapped to a different chr     #PE read中，两端分别比对到两条不同的序列的reads总数。
1712129 + 0 with mate mapped to a different chr (mapQ>=5)     #PE read中，两端分别比对到两条不同的序列，且mapQ>=5的reads总数。

根据名称从fasta文件提取序列

faidx

samtools faidx input.fasta

genome.fasta.fai

column 1 : header name
column 2 : sequence size ( how many nucleotid)
column 3 : file offset where the sequence start. ( How many caracter from the top of the file )
column 4 : line size ( how many nucleotid per line without cariage or line return )
column 5 : full line size ( how many caracter with line return )

samtools faidx genome.fasta chr1  > chr1.fasta

查看alignments数量(注意不是reads的数量)

samtools view -c test.bam

$ samtools view  hisat2.bam | wc -l
132342474
$ samtools view  -c hisat2.bam 
132342474

depth

统计一号染色体上的比对的各位点覆盖深度

samtools depth -r chr1 test.bam

bam to fastq

samtools \
    bam2fq \
    -@ 2 \
    test2-test_L1.bam | gzip --no-name > test2-test_L1_interleaved.fq.gz

samtools fastq \
    -1 fastq/reads1.fastq \
    -2 fastq/reads2.fastq \
    -s fastq/singletons.fastq \
    2-germline/bams/mother.bam

• 其中fastq/reads1.fastq与fastq/reads2.fastq是输出的reads1和reads2文件的名称；
• fastq/singletons.fastq是包含未匹配的single-end reads的输出文件的名称
• 2-germline/bams/mother.bam是输入的BAM文件的名称

从fasta文件提取指定区域的序列

samtools faidx output/index/ref.fasta  20:1000000-1000200

提取指定区域或染色体上比对的reads

samtools index   hisat2.bam
samtools view   -b  hisat2.bam   12 | samtools fastq  |  gzip  >  gwat_ABX_OLZ_18.fq.gz
# samtools view   -b  hisat2.bam   12 | samtools bam2fq   |  gzip  >  gwat_ABX_OLZ_18.fq.gz

samtools index   hisat2.bam
samtools view  -f 66  -b  hisat2.bam   12 | samtools fastq  | gzip  >  gwat_ABX_OLZ_18.fq_1.gz

https://www.samformat.info/sam-format-flag

获取提取到的reads的名称

zcat  gwat_ABX_OLZ_18.fq_1.gz |  awk 'NR%4==1 {print  substr($1,2)}'  > reads.txt
sed -i  "s/\/1/\/2/" reads.txt

接下来就是根据reads名称提取reads2，并且将reads的顺序和reads2的顺序调整为一致

samtools view  -f 130  -b  hisat2.bam    12 | samtools fastq  | seqtk  subseq -  reads.txt  | gzip  >  gwat_ABX_OLZ_18.fq_2.gz

process BAM_FASTQ{
    publishDir = [
        path: { "test_data/fastq" },
        mode: 'copy',
        saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
    ]
    debug true
    label 'process_high'

    
    input:
    tuple val(pair_id),path(bam)
    
    output:
    tuple val(pair_id),path("${pair_id}/*.gz")

    script:
        // println reads[0]
    """
        [ ! -d ${pair_id} ] && mkdir ${pair_id}
        samtools index   ${bam} 
        samtools view  -f 66  -b  ${bam}   12 | samtools fastq  | gzip  >   ${pair_id}/${pair_id}_1.gz
       zcat  ${pair_id}/${pair_id}_1.fastq.gz |  awk 'NR%4==1 {print  substr($1,2)}'  > reads.txt
        sed -i  "s/\/1/\/2/" reads.txt 
        samtools view  -f 130  -b  ${bam}    12 | samtools fastq  | seqtk  subseq -  reads.txt  | gzip  >    ${pair_id}/${pair_id}_2.gz
    """

}

workflow{
    BAM_FASTQ(["test","/ssd1/wy/workspace/RNA-seq/workspace/results/hisat2/gwat_ABX_OLZ_18/hisat2.bam"])
}

mpileup

mpileup是samtools的一个命令，用来生存bcf文件，然后再用bcftools进行SNP和Indel的分析。另外，bcftools是samtools的附带软件。

https://www.jianshu.com/p/1878ebd9930a

对bam文件进行排序

samtools sort EV_str.bam -o EV_strain.sort.bam

一次分析所有的数据

hisat2 -x ../../data/reference/yeast_ref -U ../../data/raw_data/SRR1916152.fastq | samtools view -bS - | samtools sort - -o EV_3.bam
hisat2 -x ../../data/reference/yeast_ref -U ../../data/raw_data/SRR1916153.fastq | samtools view -bS - | samtools sort - -o EV_4.bam
hisat2 -x ../../data/reference/yeast_ref -U ../../data/raw_data/SRR1916154.fastq | samtools view -bS - | samtools sort - -o DNMT3B_2.bam
hisat2 -x ../../data/reference/yeast_ref -U ../../data/raw_data/SRR1916155.fastq | samtools view -bS - | samtools sort - -o DNMT3B_3.bam
hisat2 -x ../../data/reference/yeast_ref -U ../../data/raw_data/SRR1916156.fastq | samtools view -bS - | samtools sort - -o DNMT3B_4.bam

使用samtools tview查看比对信息

samtools tview EV_3.bam

depth可以得到每一个碱基位点上的测序深度

samtools depth samp.bam > base_depth.txt
sambamba depth -t 10 -w 500 samp.bam > 500base_depth.txt
# -w breadth of the window, in bp, 计算窗口的平均深度

bedtools  genomecov -ibam input.bam -d > sample.depth

统计某个染色体上的某个区域的平均深度

chr=chr6 && start=52156548 && end=52161782 && samtools depth -r $chr:$start-$end test.bam | awk  -v start=$start -v end=$end '{sum+=$3} END { print "Average = ",sum/(end-start)}'

git clone https://github.com/shiquan/bamdst.git

需要制作一个bed文件

# tab分隔
cat >test.bed
chr6	52156548	52161782
然后

bamdst -p test.bed -o ./ test.bam
# -o：output dir

mm10.chrom.sizes

# 统计某个染色体的平均深度
grep -vi 'random' mm10.chrom.sizes  | grep -vi un | while read i;do 
	conf=($i)
	chr=${conf[0]}
	len=${conf[1]}
	samtools depth -r $chr test.bam | awk -v chr=$chr -v len=$len '{sum+=$3} END { print chr,"-Average = ",sum/len}';done
# chr1 -Average =  1.01466
# chr2 -Average =  1.06905

samtools idxstats --threads 10  workspace/results/hisat2/control1/hisat2.bam
1       195154279       7006995 494986
2       181755017       5763090 524057
3       159745316       4296878 491237
4       156860686       7296760 628701
5       151758149       5953368 412401

其中第一列是染色体名称，第二列是序列长度，第三列是mapped reads数，第四列是unmapped reads数。

技巧五、质量值为30、测序深度大于10的一号染色体上的比对信息

samtools view -q 30 -@ 10  test.bam chr1  |samtools index|samtools depth -r chr1 |awk '$3>10{print $0}' ->q30depth10chr1.depth

flag

https://blog.csdn.net/weixin_40640700/article/details/119955063
https://www.samformat.info/sam-format-flag

• 1 ： 代表这个序列采用的是PE双端测序
• 2： 代表这个序列和参考序列完全匹配，没有插入缺失
• 4： 代表这个序列没有mapping到参考序列上
• 8： 代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
• 16：代表这个序列比对到参考序列的负链上
• 32 ：代表这个序列对应的另一端序列比对到参考序列的负链上
• 64 ： 代表这个序列是R1端序列， read1;
• 128 : 代表这个序列是R2端序列，read2；
• 256： 代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的
• 512： 代表这个序列在QC时失败了，被过滤不掉了（# 这个标签不常用）
• 1024: 代表这个序列是PCR重复序列（#这个标签不常用）
• 20q48: 代表这个序列是补充的比对（#这个标签具体什么意思，没搞清楚，但是不常用）

samtools flags 4
# 0x4     4       UNMAP

查找flag包含4(序列没有mapping到参考序列上)的比对结果

samtools view  -f 4  hisat2.bam | less -S 

E100050738L1C014R02202463985    133     1       9503    0       *       =       9503    0       CAATAAAGATTTCACCTTTAATAATACTAGGTGGAATCTTTATTGGAGTTTGGTTATCAACTTTCCAAGAAAGGGTCATGTATTTGAAGATTGAACCCTGGCATGTGGCCCTGTTTGGGCAAGTTGTAGGAAATCTGGGAAATATAGCTT FFEFFFCFFFFFDFFCFDFFFFFFEFFFFFEEFFFEFEFFFFFFFEEEEEFEADFFEFDFAFFFEFFFFEFE<<FD>FFFFFFFFFFFEEEFFEEEBFEFBFFEFEGEFBED>FFFFFDFD=DDFFEFFEEF=D@FFBECEFE>ECD@DF YT:Z:UP
E100050738L1C036R04102024167    69      1       10633   0       *       =       10633   0       ATCGTTCAAAGTATTTTTTTCTCTATCTTTTCAAATATATTAATCCACTGCGTTTTGTGATGCTAGCTTTGTTTTTTTTTTTTTTTTTTTTTGGTTTTTTTTTTTGGAGGTGGTGGTTTTGATTGGTGGTATGGCGCTGTCCTACGAGAT EE>GEDFEFFFE?FEFF1F?FFFF2FDFEFC?FDFF:FDFDDDE@DDDF>2(DF0F@6>B,FCE>;*E3E+EEDFEE6ECF9EFCBFFF@FE/B.)*E:E?AF:5*A/9/;.;=C(,EC7C+0)+>/)F2>1*(*73E:3?BDD<*80/4 YT:Z:UP

使用samtools flag 133查看flag值为133的含义是PAIRED,UNMAP,READ2

使用samtools view -F 4 hisat2.bam查找flag值不包括4的比对结果

https://zhuanlan.zhihu.com/p/545946218
https://broadinstitute.github.io/picard/explain-flags.html
https://www.jianshu.com/p/ff6187c97155