bam文件转fq.gz文件

S_AGZX

已于 2022-06-19 19:06:53 修改

阅读量2.9k

点赞数

分类专栏：测序基础文章标签： anaconda

于 2022-06-19 19:06:16 首次发布

本文链接：https://blog.csdn.net/s_agzx/article/details/125312146

版权

测序基础专栏收录该内容

5 篇文章 0 订阅

订阅专栏

bam转fq过程

1、fq文件格式
2、bam文件格式
3、转换思路
4、参考资料

1、fq文件格式

fastq格式是一种包含质量值的序列文件，一般用来存储原始测序数据。下面是fastq格式常见的序列格式。

@FCD056DACXX:3:1101:2163:1959#TCGCCGTG/1
TCCGATAACGCTCAACCAGAGGGCTGCCAGCTCCGATCGGCAGTTGCAACCCATTGGCCGTCTGAGCCAGCAACCCCGGA
+
gggiiiiiiiiiiiiiiiiiiiiiiiiiigggggeeecccccc^bcbcccccccbccccc]aaccbbccc^R^^acccc_
@FCD056DACXX:3:1101:2194:1984#TCGCCGTG/1
AGACGACGACTTCGTTTCCCGCCGCGAGTTGCGCCATGATCGCGGTGTGCAGATTCGTTACGCCCTGGGCCACGGAGACG
+
gggiihiiiiiiiihiiiiiiiiiigeccccccccccccccccccaccccdcccccccccccacc_accccccccccV^^

第一行：以‘@’开头，是这一条read的名字，这个字符串是根据测序时的状态信息转换过来的，中间不会有空格，它是每一条read的唯一标识符，同一份FASTQ文件中不会重复出现，甚至不同的FASTQ文件里也不会有重复；
第二行：测序read的序列，由A，C，G，T和N这五种字母构成，这也是我们真正关心的DNA序列，N代表的是测序时那些无法被识别出来的碱基；
第三行：以‘+’开头，在旧版的FASTQ文件中会直接重复第一行的信息，但现在一般什么也不加（节省存储空间）；
第四行：测序read的质量值，这个和第二行的碱基信息一样重要，它描述的是每个测序碱基的可靠程度，用ASCII码表示。

2、bam文件格式

bam、sam文件内容其实是一样的，只是bam是二进制的压缩文件，需要通过特定的软件来进行查看。
bam文件通常分为header section（头部分，注释信息，以@开头，可有可无）和alignment section（比对结果）两个部分。比对结果由12个字段组成。
示例如下：

F300000235L1C002R0530799573_TGCTTCCTACCA    1187    chr13   32968914    60  82M9D11M    =   32969045    224 GCCTCATATGTTAATTGCTGCAAGCAACCTCCAGTGGCGACCAGAATCCAAATCAGGCCTTCTTACTTTATTTGCTGGAGATTTTTCTGCTAG   A/B7>FDCFDEFDBFEBDFCAEFEAAACAECEDE@B@=B?5BDF1D;:EFCE5ADE>BEE;CCDF=FF=FEBEBEEE=EE@D@AC=<DCACDB   X0:i:1  X1:i:0  MD:Z:82^TTTTCTGTG11 PG:Z:MarkDuplicates RG:Z:case   XG:i:9  DI:i:6783411    AM:i:37 NM:i:9  SM:i:37 XM:i:2  XO:i:1  MQ:i:60 DS:i:2  XT:A:U

其中：
（1）F3*_TGCTTCCTACCA：序列的名字，也就是reads的名称
（2）118： FLAG值，是一个标记的数字，均可以转换为2的n次方。这种特定能保证它们数字加和是唯一的。
可以使用 samtools flags 1147 了解flag值由哪些数字组成，或者到网站了解详细含义：http://broadinstitute.github.io/picard/explain-flags.html。

1：序列是一对序列中的一个；
2：比对结果是一个pair-end比对的末端；
4：没有找到位点；
8：这个序列是pair中的一个但是没有找到位点；
16：在这个比对上的位点，序列与参考序列反向互补；
32：这个序列在pair-end中的的mate序列与参考序列反响互补
64：序列是 mate 1；
128：序列是 mate 2；

假如说标记不是以上列举出的数字，比如说83=（64+16+2+1），其实就是这几种情况值和，由二进制混合来表示。

（3） chr13：参考序列的名字
（4）32968914：在参考序列上的位置
（5） 60：mapping qulity 越高则位点越独特，比对的质量值。
比对的质量值 0-60 依次表示未比对上-最高比对值，255表示结果不可用。
（6）85M9D8M：CIGAR 值
表示reads的比对状态，M表示匹配，I插入，D缺失，N可变剪切，S表示替换，H表示切除，例如这里就是这条reads前面85个碱基比对上了，然后中间9个缺失，后面8个又比对上了；
更详细的可以见生物信息学100个基础问题 —— 第19题 SAM/BAM中的CIGAR值

Soft Clip，是指虽然比对不到基因组，但是还是存在于SEQ (segment SEQuence)中的序列，此时CIGAR列对应的S(Soft)的符号。直白点说，就是虽然比对不上参考基因组，但是在BAM/SAM文件中的reads上还是存在的序列（并没有被截断扔掉的序列）

Hard Clip，同样的，就表示比对不上并且不会存在于SAM/BAM文件中的序列（被截断扔掉了的序列，此时CIGAR列会留下H(Hard)的符号，但是序列的那一列却没有对应的序列了）。

只要一条序列两端有比对不上的序列部分，就是Soft Clip，这个一条序列上有比对不上的部分的现象是必然存在的（比如结构变异的断点的部分），这种两端比对不上的read的特殊的表示方法，就是Soft Clip。Soft Clip是可以独立存在的。

Soft Clip不一定有Hard Clip，而有Hard Clip则一定有Soft Clip。Hard Clip存在的本意，是减少BAM文件序列的冗余度，比如有一条read，它能比对到A，B两个地方，在A地方，是60M90S，在B地方是60H90M，此时一条read其实已经在A位置有了完整的序列信息，在B位置的信息其实是冗余的，所以在B位置可以引入Hard Clip这样一个标记形式，就能把B位置的序列标记为secondary。

（7）=（是否存在多重比对）：mate 序列所在参考序列的名称，mate一般指大的片段序列
reads二次比对结果，=表示没有二次结果；如果3和7都是*，表示没有结果；如果只有7是*，表示只比对到3的结果，如下：

chr13	32968934	37	62M9D31M	chr12	52009363	0

（8）32969045：mate 序列在参考序列上的位置，即多重比对时另一条reads位置；
（9）224：文库插入片段长度，当mate 序列位于本序列上游时该值为负值。
（10）read的序列：

GCCTCATATGTTAATTGCTGCAAGCAACCTCCAGTGGCGACCAGAATCCAAATCAGGCCTTCTTACTTTATTTGCTGGAGATTTTTCTGCTAG

（11）read序列对应的ASCII码格式的碱基质量值：

A/B7>FDCFDEFDBFEBDFCAEFEAAACAECEDE@B@=B5BDF1D;:EFCE5ADE>BEE;CCDF=FF=FEBEBEEE=EE@D@AC=<DCACDB

（12）可选的区域 header section：

X0:i:1  X1:i:0  MD:Z:82^TTTTCTGTG11 PG:Z:MarkDuplicates RG:Z:case   XG:i:9  DI:i:6783411    AM:i:37 NM:i:9  SM:i:37 XM:i:2  XO:i:1  MQ:i:60 DS:i:2  XT:A:U

其中header section用不同的tag表示不同的信息，主要有@HD，说明符合标准的版本、对比序列的排列顺序；@SQ，参考序列说明；@RG，比对上的序列（read）说明；@PG，使用的程序说明；@CO，任意的说明信息。Tag以键值对的形式存在。

	AS:i 匹配的得分
	XS:i 第二好的匹配的得分
	YS:i mate 序列匹配的得分
	XN:i 在参考序列上模糊碱基的个数
	XM:i 错配的个数
	XO:i gap open的个数
	XG:i gap 延伸的个数
	NM:i 经过编辑的序列
	YF:i 说明为什么这个序列被过滤的字符串
	YT:Z
	MD:Z? 代表序列和参考序列错配的字符串

3、转换思路

3.1 软件bedtools自带功能

主要是以参考资料中的信息为主，大致步骤如下：
（1）对bam文件重排序
流程跑出来的input.bam是以染色体位置排序的，重新sort成按reads名称排序的形式；如果本来就是，这步可跳过。

#sort 
samtools sort -n input.bam -o input.name.bam

（2）转换成fastq文件
使用bamtofastq功能把reads还原成fq文件。需注意，仅有bam文件中read_1、read_2都存在（两条reads排序后挨着）的reads才会被保留，仅有1条的会被舍弃。

#convert  
bedtools bamtofastq -i input.name.bam \
-fq out.R1.fq \
-fq2 out.R2.fq

（3）压缩fastq文件
gzip压缩成fq.gz形式。

#gzip fq
gzip out.R1.fq
gzip out.R2.fq

3.2 自己写代码

（1）找bam => fq转化规律
通过flag值得到两组指标READ1 or READ2, MREVERSE or REVERSE;
其中，READ1-READ2来决定放到哪个fastq文件中，MREVERSE-REVERSE看测序序列是否与参考序列相同(MREVERSE=相同)。MREVERSE时bam中信息与fq一致, 无需改变; REVERSE时需要质量值反向、序列反向互补;

（2）脚本实现逻辑
模拟后的bam(sort_by_pos) => 读入pysam => 保留flag值<2048的reads(去掉补充序列) => 按1中规则转换，同时修改reads名称格式,dict记录 => 仅保留reads有read1+read2信息的那种 => 写入并压缩fq文件;

3.3 代码示例

# !D:\00_software\13.Python\python.exe
# -*- coding: UTF-8 -*-
# @File    : bam_to_fq.py

import argparse
import os
import sys
import pysam
from collections import defaultdict


def get_opt():
    group = argparse.ArgumentParser(description="this script is convert file format. bam(sort_by position) => fq.gz")
    group.add_argument("-i", "--input_file", help="input bam file", required=True)
    group.add_argument("-n", "--file_name", help="specify the result file(fq.gz) name", required=True)
    group.add_argument("-o", "--out_dir", help="output directory", default="./")

    return group.parse_args()


def sort_bam(bam_by_position):
    samtools = "/data/Apps/Production/samtools-1.7/bin/samtools"
    is_exists = os.path.exists(samtools)
    if not is_exists:
        print("%s 's abs_path is invalid, please check!!!" % samtools)
        sys.exit(1)

    bai_file = bam_by_position + ".bai"
    if not os.path.exists(bai_file):
        cmd_bai = samtools + " index " + bam_by_position
        os.system(cmd_bai)

    prefix_bam = bai_file.split(".bam")[0]
    sort_bam_file = prefix_bam + "_sort.bam"
    if os.path.exists(sort_bam_file):
        os.system("rm {}".format(sort_bam_file))
    temp_path_list = os.popen("ls {}*sort.bam.tmp*".format(prefix_bam))
    for path_list in temp_path_list:
        os.system("rm {0}".format(path_list))

    cmd_sort = samtools + " sort -n " + bam_by_position + " -o " + sort_bam_file
    os.system(cmd_sort)
    return sort_bam_file


def str_complement(sequence):
    # 构建互补字典
    comp_dict = {
        "A": "T",
        "T": "A",
        "G": "C",
        "C": "G",
        "a": "t",
        "t": "a",
        "g": "c",
        "c": "g",
        "N": "N"
    }
    #求互补序列
    sequence_list = list(sequence)
    sequence_list = [comp_dict[base] for base in sequence_list]
    string = ''.join(sequence_list)
    return string


def str_reverse(sequence):
    return sequence[::-1]  # 求反向序列


def reads_tag(input_file):
    samfile = pysam.AlignmentFile(input_file, "rb")
    # samfile = pysam.AlignmentFile("test.bam", "rb")
    all_reads_dict = defaultdict(dict)
    for read in samfile.fetch():
        # 过滤掉补充reads
        if read.is_supplementary:
            continue

        if read.is_read1:
            read_orientation = "read1"
            read_name = "@" + read.qname + " 1:N:0:AGCGATAG+ACGAATAA"
        else:
            read_orientation = "read2"
            read_name = "@" + read.qname + " 2:N:0:AGCGATAG+ACGAATAA"
        fixed = "+"
        raw_read_name = read.qname

        if read.is_reverse:
            temp_seq_1 = read.seq
            temp_seq_2 = str_reverse(temp_seq_1)
            read_seq = str_complement(temp_seq_2)

            temp_quality = read.qual
            read_quality = str_reverse(temp_quality)
        else:
            read_seq = read.seq
            read_quality = read.qual
        total_reads = "{name}\n{seq}\n{fix}\n{qual}\n".format(name=read_name, seq=read_seq, fix=fixed, qual=read_quality)
        # print(read.qname, read.is_reverse, read.seq, read.qual)
        # print("total_reads", total_reads)
        all_reads_dict[raw_read_name][read_orientation] = total_reads
    return all_reads_dict


def save_fq(all_reads_dict, file_name, out_dir):
    fq1 = out_dir + "/" + file_name + "_R1.fq"
    fq2 = out_dir + "/" + file_name + "_R2.fq"
    ff_1 = open(fq1, 'w')
    ff_2 = open(fq2, 'w')

    for read_name in all_reads_dict.keys():
        orient_num = len(all_reads_dict[read_name].keys())
        # 只保留read_name有两条信息的测序数据
        if orient_num != 2:
            continue
        else:
            # print(read_name, orient_num)
            read1_info = all_reads_dict[read_name]['read1']
            read2_info = all_reads_dict[read_name]['read2']
            ff_1.write(read1_info)
            ff_2.write(read2_info)
    ff_1.close()
    ff_2.close()

    fq_gz_1 = out_dir + "/" + file_name + "_R1.fq.gz"
    fq_gz_2 = out_dir + "/" + file_name + "_R2.fq.gz"
    if os.path.exists(fq_gz_1) or os.path.exists(fq_gz_2):
        os.system("rm {} {}".format(fq_gz_1, fq_gz_2))

    os.system("/usr/bin/gzip {}".format(fq1))
    os.system("/usr/bin/gzip {}".format(fq2))
    return fq_gz_1, fq_gz_2


def main():
    opts = get_opt()
    input_file = opts.input_file
    file_name = opts.file_name
    out_dir = opts.out_dir
    # print(input_file, file_name, out_dir)

    all_reads_dict = reads_tag(input_file)
    fq_gz_1, fq_gz_2 = save_fq(all_reads_dict, file_name, out_dir)
    print(fq_gz_1, fq_gz_2)


if __name__ == '__main__':
    main()