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:ZMD: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.pyimport argparse
import os
import sys
import pysam
from collections import defaultdictdef 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_positionos.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_fileos.system(cmd_sort)return sort_bam_filedef 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 stringdef 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():# 过滤掉补充readsif read.is_supplementary:continueif 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.qnameif read.is_reverse:temp_seq_1 = read.seqtemp_seq_2 = str_reverse(temp_seq_1)read_seq = str_complement(temp_seq_2)temp_quality = read.qualread_quality = str_reverse(temp_quality)else:read_seq = read.seqread_quality = read.qualtotal_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_readsreturn all_reads_dictdef 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:continueelse:# 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_2def main():opts = get_opt()input_file = opts.input_filefile_name = opts.file_nameout_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()

4、参考资料

（1）bamtofastq 将bam文件转成fastq文件
（2）bam格式转换为Fastq/Fasta格式

bam文件转fq.gz文件相关推荐

Linux读取下机数据.fq.gz文件
#读取压缩文件 $ zcat R1.fq.gz 会唰唰唰全部读取,按Ctrl + C终止 #只查看部分 $ less R1.fq.gz 按↓键可往下翻页,按q退出
linux上怎么解压zip文件和tar.gz文件
解压zip文件的方法使用命令: unzip xxx.zip 解压tar.gz文件的方法使用命令: tar -zxvf xxx.tar.gz
python读取nii文件、nii.gz文件
显示标准nii.gz或nii文件 import numpy as np import nibabel as nib from ipywidgets import interact, interacti ...
windows环境下安装Python的.whl文件和tar.gz文件
一.whl文件的安装: 1.先弄清楚自己的Python是什么版本的,以方便后续下载合适的.whl文件. win+R进入命令运行窗口,输入cmd打开命令提示符,接着输入python即可 2.选择需要的 ...
.zip 文件和 .tar.gz文件的区别
经常去网站下载资源,看到后缀名,很疑惑,遂,今查之. 简单来说: tar.gz压缩格式用于unix的操作系统,而zip用于windows的操作系统,但在windows系统中用WinRar工具同样可以解 ...
Nibabel 读取 nii 文件和 nii.gz 文件
读取nii文件并且将nii文件转换为png格式 import numpy as np import nibabel as nib import os import imageio# 文件路径 nii_ ...
CT图像分割dicom文件与nii.gz文件预处理----窗宽（window width）和窗位（window level）的设置
最近被CT图像的值弄得很烦,记录一下. CT分割也是个很热门的话题,病灶分割,器官分割等. CT图像大多是两种格式.dcm和nii.gz,当然也有别的,但这里我就不说别的,就说这两种常用的. .dcm ...
解压tar.xz文件和tar.gz文件
tar xvJf *tar.xz tar zxvf *tar.gz
linux 解压.tar.gz文件
(1)解压 .tar.gz文件 tar -zxvf 文件名.tar.gz 其中,文件名.tar.gz 是你要解压的文件的名称. 解释一下命令的选项: -z:表示使用 gzip 压缩算法进行解压. -x ...

bam文件转fq.gz文件

bam转fq过程