肿瘤生信数据分析知识积累

bed 文件

probe bed 和 target bed:

bed 文件是什么？

bed文件是记录基因组位置信息的标准文件格式，同时也用于存储与位置相关的信息，例如在ChIP-Seq 分析中，长以bed文件存储检测信号强度的信息、结构变异检测（SV）结果也可以用bed文件或bedpe文件进行存储。

bed 文件格式

第一列：位置所在的染色体
第二列： 表示区域在染色体的起始位置
第三列：表示区域在染色体上的终止位置
第四列： 基因名/注释信息
第五列： 比对的参考基因组是正链还是负链。
还有些其他的可选信息，未列出。

fastq 文件 测序数据的质量信息

FASTQ是基于文本的，保存生物序列（通常是核酸序列）和其测序质量信息的标准格式。其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。

• 格式说明
  FASTQ文件中每个序列通常有四行：

第一行：序列标识以及相关的描述信息，以‘@’开头；
第二行是序列
第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加
第四行，是质量信息，和第二行的序列相对应，每一个序列都有一个质量评分，根据评分体系的不同，每个字符的含义表示的数字也不相同。

fastqc 软件

FastQC是一款基于Java的软件，它可以快速地对测序数据进行质量评估，其官网为：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

我们可以输入 fastqc -h来查看fastqc所有的命令和参数

一般最常用的是
fastqc file1 -o file2

fastqc 报告解读

可以去官网看。
绿色代表合格，黄色代表警告，红色代表不合格。

filename即为被质检的文件名
file type指文件类型
encoding指测序平台的版本和相应的编码版本号
total sequences指reads的数量
sequence length指测序长度
%GC表示整体序列的GC含量；二代测序深度高，GC含量偏高。

per base sequence quality显示的是每个碱基的质量
横轴是各个碱基的位置，我这里只有第一个到第150个，纵轴是质量得分
Q = -10*log10（error P）
即当值为30时，其Q等于0.001，表示0.001的错误率，所以当纵轴的值越大时，其质量越好。一般都在30以上，说明该条序列质量良好。可以看到这个例子后面的测序质量较好。

该图显示的是每条read的平均质量
横轴为上副图的纵轴，纵轴为read的数量，可以看到该条序列大部分read的质量都在30左右，所以质量还可以。

该图显示的是各个碱基的含量。
横轴代表各个碱基的位置，纵轴表示百分比
A和T应该相等，C跟G应该相等，但一开始不稳定，所以可能会出现开头的情况。

这幅图代表GC含量，横轴代表GC含量百分比，纵轴表示read数量
蓝色为理论值，红色为实测值，两者应越接近越好。

改图表示的为reads的N含量。
当fastqc无法识别某个位置的碱基为ATGC时，就会产生N。

改图显示的为测序长度
横轴代表长度，纵轴代表数量。理论上每条长度应该一致，但难免会出现偏差。该例子中大部分长度都为150。

改图显示的是完全相同序列的重复次数
横轴代表重复的次数，纵轴代表百分比含量

最后一张图表示的每个碱基位置上的adapter含量。

Q20 Q30 计算

软件 Q30

碱基质量分数和错误率是衡量测序质量的重要指标，质量值越高代表碱基被测错的概率越小。
行业中Q20与Q30则表示质量值≧20或30的碱基所占百分比。例如一共测了1G的数据量，其中有0.9G的碱基质量值大于或等于20，那么Q20则为90%。

Q30代表碱基的正确判别率是99.9%，错误率为0.1%。同时我们也可以理解为1000个碱基里有1个碱基是错误的。

Q20代表改位点碱基的正确率是99%，错误率为1%。

对于整个数据来说，我们认为每100个碱基里可能有一个是错误的，在碱基质量模块报告的坐标图里，背景色沿着y-轴将坐标图氛围3个部分；最上面的绿色是碱基质量很好的区，Q值在30以上。中间的橘色是碱基质量在一些分析中可以接受的区，Q值在20-30之间。最下面的红色是碱基质量很差的区。
检查差异表达为目的的RNA-seq分析，一般要求碱基质量在Q20以上就可以了。
但在以检查变异为目的的数据分析中，一般要求碱基质量要在Q30以上。

测序质量分数的分布的两个特点：
1.测序质量分数会随着测序循环的进行而降低。
2.有时每条序列前几个碱基的位置测序错误率较高，质量值相对较低。

质量值是Q20，则错误识别的概率是1%，即错误率1%，或者正确率是99%；
质量值是Q30，则错误识别的概率是0.1%，即错误率0.1%，或者正确率是99.9%；
质量值是Q40，则错误识别的概率是0.01%，即错误率0.01%，或者正确率是99.99%；

Q20Q30报告：
Total Bases: 总碱基数
Q20 Bases:Q20的碱基数
Q30 Bases:

fastp 软件进行 质控和过滤

执行FASTQ数据的质量控制，支持单端和双端短读数据，作为一体化的FASTQ预处理器，fastp提供了诸如质量分析，接头修整，读取过滤和基本校正等功能。在处理完所有读取后，这些值将合并，并且报告者将生成HTML和JSON格式的报告。
$fastp -q 15 -u 40 -3 -W 4 -M 20 -n 5 -l 15 -w $nt -i $R1 -I $R2 -o ${sample}.R1.paired.fastq.gz -O ${sample}.R2.paired.fastq.gz

这是一个Shell脚本中的命令，使用了fastp工具对输入的两个fastq文件进行质量控制和过滤，并将过滤后的结果保存到$sample.R1.paired.fastq.gz和${sample}.R2.paired.fastq.gz中。具体来说，该命令使用了以下参数：
-q 15：设定质量阈值为15，低于该阈值的reads将被过滤掉。
-u 40：设定5’端截断长度为40，即将每个read的前40个碱基去除。
-3：对3’端进行质量截断，默认截断位置为质量值低于15的第一个碱基。
-W 4：设定滑动窗口大小为4，用于计算平均质量。
-M 20：设定最小平均质量阈值为20，低于该阈值的reads将被过滤掉。
-n 5：设定最多允许5个N碱基存在于read中，超过该数量的reads将被过滤掉。
-l 15：设定最短read长度为15，短于该长度的reads将被过滤掉。
-w $nt：设定线程数为
n t，即使用nt，即使用nt个线程并行处理。
-i $R1：指定输入的R1文件。
-I $R2：指定输入的R2文件。
-o ${sample}.R1.paired.fastq.gz：指定输出的R1 paired reads文件。
-O ${sample}.R2.paired.fastq.gz：指定输出的R2 paired reads文件。

该图是${sample}.R1.paired.fastq.gz 的文件内容。

数据比对 bwa

比对的目的：测序仪的原理，边合成变测序，测序时，先把RNA逆转录成双链的cDNA, 有一条链与原来的模版是相同的，另一条链与原来的模板反向互补的。普通转录本测序的时候是对两条链都进行测序，打断后得到些reads，这些reads 有的来自模板链，有的来自反向互补链，对所有的reads都进行测序的叫普通转录本测序。

链特异性测序： 给要测序的链进行加标记，这样测序的时候，就能把有特异性标记的reads找出来。就是只比对+链的reads。

比对，就是把测序的小片段reads与参考基因组进行比较，大概推测我们的这个基因来自哪里。

双端测序比对：

单端测序比对：

BWA是一款基于BWT的快速比对工具，其由三个算法组成。这三个算法分别是：BWA backtrack, BWA SW and BWA MEM。其中，BWA MEM是最新的，其更快更准确，更适合用于人重数据分析。对于上述三种算法，首先需要使用索引命令构建参考基因组的索引，用于后面的比对。所以，使用BWA整个比对过程主要分为两步，第一步建索引，第二步使用BWA MEM进行比对。

对于参考基因组文件大于2G的使用bwtsw算法，使用bwtsw算法必须保证参考基因组文件大小大于10M。
$bwa mem -M -R "@RG\tID:${sample}\tLB:${sample}\tSM:${sample}\tPL:$pl" -t $nt $hg19 ${sample}.R1.paired.fastq.gz ${sample}.R2.paired.fastq.gz > ${sample}.sam

    参数详解：R——设置reads标头，“\t”分割。 例如：’@RG\tID:foo\tSM:bar’；

                       t——设置线程数；

                       M——将较短的split hits标记为secondary，与picard兼容。
这个例子中：
${sample}： 一个文件夹的名字—— N19-01867E1L1_raw
$nt ： 线程数 - 10
$hg19： hg19的路径- 人类基因组
${sample}.R1.paired.fastq.gz： 要比对的数据1
${sample}.R2.paired.fastq.gz： 要比对的数据2
${sample}.sam： 比对好的数据保存的地址

${sample}.R1.paired.fastq.gz 代表： N19-01867E1L1_raw.R1.paired.fastq.gz
${sample}.R1.paired.fastq.gz 的内容是：

/home/biosoft/bwa-0.7.17/bwa mem -M -R @RG\tID:N19-01867E1L1_raw\tLB:N19-01867E1L1_raw\tSM:N19-01867E1L1_raw\tPL:ILLUMINA -t 10 /data/database/1000g_genome/hg19_1000g/hg19_1000g.fa N19-01867E1L1_raw.R1.paired.fastq.gz N19-01867E1L1_raw.R2.paired.fastq.gz

sam文件格式解析：
官方解读： http://samtools.github.io/hts-specs/SAMv1.pdf

SAM格式文件一般是序列比对程序标准输出文件，以TAB作为分隔符，并且前几行通常是一些header（可选）

1. 文件中的header行
   SAM文件中的header以 @开头
   每个SAM文件的开头可以使用 @HD VN用于指定当前sam文件的版本信息
   例如@HD VN:1.0 SO:unsorted就表示当前SAM文件格式版本为1.0，SO表示SAM文件是否为排序好的

比对结果的可视化：

比对质量表示： MAPQ，即比对的质量，一致性如何。

测序深度和覆盖度：
https://www.bilibili.com/video/BV1GM411a74W/?spm_id_from=333.337.search-card.all.click&vd_source=b9a2c87aaaf711b6e400c3041097d8ff


因为只有绿色的部分被覆盖到了，所以参考基因组的覆盖度就是46/112=41.07%

Fig.2 比对部分解释

（1）第一列：QNAME

reads序列的名称，如果QNAME唯一，则序列被认为来源于同一模板；‘*’表示该字段缺省；一般情况下，该字段为FASTQ文件的第一行信息；嵌合（Chimeric alignment）比对或者多次比对（Multiple mapping）的序列会导致一个QNAME在SAM中多次出现。

（2）第二列：FLAG

二进制标志位，表示读取序列的比对状态和其他信息。可使用第三方网站输入数值查询具体含义。

Fig.3 FLAG质量值含义

（3）第三列：RNAME

染色体名称。如果这列是“ * ”，可以认为这条read没有比对上的序列，则这一行的第四，五，八，九 列是“0”，第六，七列与该列是相同的表示方法。

（4）第四列：POS

比对的位置，从对应上的染色体第1位开始往后计算。该read如果没有比对上参考序列，此处为0且RNAME和CIGAR也无值。

（5）第五列：MAPQ

MAPQ比对质量值，值越高说明该read比对到参考基因组上的位置越唯一，计算公式为：-10logP{错比概率} 。一般结果是这一列的数值是从0到60，且0和60这两个数字出现次数最多，其中0表示在参考基因组有多种定位的可能性，60表示在参考基因组只有这一种定位位置。

（6）第六列：CIGAR

表示比对过程中的操作，如匹配、插入或删除等。比如，一条150bp长的read比对到基因组之后，假如看到它的CIGAR字符串为：33S117M，其意思是说在比对的时候这条read开头的33bp在被跳过了（S），紧接其后的117bp则比对上了参考序列（M）。这里的S代表软跳过（Soft clip），M代表匹配（Match）。N表示可变剪接位置，常见于RNA-seq。H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现。

Fig.4 标记字符含义

（7）第七列：RNEXT

双端测序中另外一条read比对的参考序列的名称，单端测序此处为0。=表示参考序列与reads一模一样，*表示没有完全一模一样的参考序列。

（8）第八列：PNEXT

双端测序中，是指另外一条read比对到参考基因组的位置坐标，最小为1（1-based leftmost）

（9）第九列: TLEN

序列模板长度，如果同一个片段都比对上了同一个参考序列，为最左边的碱基位置到最右边的碱基位置（左为正，右为负）。当mate 序列位于本序列上游时该值为负值。不可用时，为0。

（10）第十列：SEQ

表示测序得到的read序列, FASTQ的第二行。

（11）第十一列：QUAL

read序列的质量值，表示每个碱基的测序质量值,FASTQ的第四行。

（12）第十二列：Optional fields

可选的区域。格式类似AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:35T0 YT:Z:UU。

AS:i 匹配的得分

XS:i 第二好的匹配的得分

YS:i mate 序列匹配的得分

XN:i 在参考序列上模糊碱基的个数

XM:i 错配的个数

XO:i gap open的个数

XG:i gap 延伸的个数

NM:i 经过编辑的序列

YF:i 说明为什么这个序列被过滤的字符串

YT:Z 值为UU表示不是pair中一部分(单末端？)、CP(是pair且可以完美匹配)

DP(是pair但不能很好的匹配)、UP(是pair但是无法比对到参考序列上)

MD:Z 代表序列和参考序列错配的字符串

GC偏好性

GC偏好性 （GC bias）

含义：在二代测序仪上测出的数据，通常都会表现出测序深度与GC 含量的相关性，称为GC bias。基因组上GC含量在50%左右的区域更容易被测到，产生的reads更多，这些区域的覆盖度更高，在高GC或者低GC区域，不容易被测到，产生较少的reads，这些区域的覆盖度更少。

备注：
1）在检测拷贝数的时候，GC含量低或者高的区域，其覆盖度小于GC含量中等的，但不意味着仅仅根据测序的覆盖度，就认为GC含量中等的拷贝数比高/低GC含量区域的高。
2）在做RNA测序分析的时候，GC含量高/低的区域reads数少，并不一定说明这个基因的表达量低。
3）在做基因组拼接的时候，因为GC偏好的存在，高/低GC含量的区域被测的少，这些区域的拼接难度就较大