人类血液 blood

人体血液主要有血细胞、血小板和血浆组成。

血细胞 Blood cell 血细胞主要包括红细胞(Red blood cell)、白细胞(White blood cell)。
血浆 Plasma 血浆(plasma)占据总血液(blood)容量的~55%，是血液中防止凝结的液体部分，更能反映血液在体内的循环情况。血浆主要是水和溶解在其中的蛋白质、激素、矿物质和二氧化碳。功能：循环液体，保持机体适宜的内环境。
血清 Serum 血清和血浆都是血液的重要组成部分，血液由血浆、血清、血小板、白细胞和红细胞组成。血浆和血清的主要区别在于它们的凝血因子不同。

1 2	血浆 = 血清 + 纤维蛋白原 + 凝血因子 Plasma = Serum + Fibrinogen + Clotting factors

血小板 Platelet 血小板是血液的一个组成部分，其功能是对血管损伤出血作出反应，形成凝块，从而形成血块。血小板没有细胞核:它们是来自骨髓巨核细胞的细胞质碎片，然后进入循环系统。
人外周血单核细胞 Human Peripheral Blood Mononuclear Cells

人外周血单核细胞(PBMC)是具有单个圆形细胞核的血细胞。这些细胞包括淋巴细胞(T、B和NK细胞)、单核细胞和树突状细胞。pbmc是免疫系统的一部分，对细胞介导和体液免疫至关重要。我们的人pbmc是从正常健康供体白细胞中分离出来的，用酸-柠檬酸-葡萄糖配方A (ACDA)作为抗凝剂收集。所有献血者必须在HBV、HCV和HIV检测中呈阴性，并得到IRB的同意。

血细胞发生(Hematopoiesis)

俗称造血，有造血干细胞完成。造血干细胞(hemopoietic stem cell)是生成各种血细胞的原始细胞,又称为多功能干细胞(multipotential stem cell)。

分离方式

采用离心分离。例如:

1.两次离心: 800Xg + 16,000Xg 用于血浆纯化(Shulin et al., Gut, 2019)。
2.两次离心: 1,600Xg + 13,000Xg, 用于一个队列收集; 2,500rpm (978Xg), 用于另外的队列收集。(Mira et al., Nature, 2022)

胞外RNA(exRNA)

细胞外的RNA(exRNA, extra-celluar RNA)包括长链和短链RNA，主要是来至血浆和血清中的RNA(cf-RNA, cell-free RNA)和富集在外泌体(exosomes)、血浆囊泡(extracellular vesicles)和血清的RNA(exoRNA)。

long RNA (> 200 nt)
mRNA, lncRNA, rRNA
small nocoding RNA (ncRNA 20-30 nt)
miRNA, piRNA, siRNA
other nocoding RNA
tRNA, Y RNA, snRNA, snoRNA, srp RNA等

Lander-Waterman 公式

Lander-Waterman formula是由Eric Lander 和 Michael Waterman于1980年提出，用于基因组组装中参数估计。

# The Lander/Waterman formula用来计算覆盖度。
# 公式的一般形式:
    C = LN / G
# C 覆盖度(coverage)
# G 单倍型基因组长度(haploid genome length)
# L 读段长度(read length)
# N 读段总数(number of reads)

参考

computational biology in Brown university Micheal Waterman web UCSD shortgun Illumina coverage

Z 算法

对于长度为n的字符串s，定义函数z[i]表示字符串本身s和其后缀s[i, n-1]的最长公共前缀长度(LCP, longest common prefix)。 z被称为s的z函数。z[0]=0。计算z数组的方法一般被称为Z-Algorithm，也称为扩增KMP算法。

Z函数的朴素实现方式朴素实现方式的复杂度为O(n^2)， for循环和while循环两层，举一个最极端的例子，s='aaaaa'，很明显每次while循环在每一个i时都会遍历全部的z，长度也为n，故朴素算法的时间复杂度为O(n^2)。

def z_function_trivial(s):
    n = len(s)
    z = [0] * n # 最长公共前缀的长度(s和其后缀)
    for i in range(1, n):
        # 索引正常 同时 *s前缀*和*s后缀的前缀*相同
        while i + z[i] < n and s[z[i]] == s[i+z[i]]:
            # 在i位置，s和s[i, n-1]的最长公共前缀的长度，相同+1
            z[i] += 1
    return z

# 举例子
# s = 'abaab'
# n = 5
#                                                           z[0] = 0
# i = 1, s = 'abaab', s的后缀s[1:4] = 'baab', 两者的LCP为0, z[1] = 0
# i = 2, s = 'abaab', s的后缀s[2:4] = 'aab' , 两者的LCP为1, z[2] = 1
# i = 3, s = 'abaab', s的后缀s[3:4] = 'ab'  , 两者的LCP为2, z[3] = 2
# i = 4, s = 'abaab', s的后缀s[4] = b       , 两者的LCP为3, z[4] = 0

Z函数线性算法

对于字符串匹配，要想提速匹配过程，必须利用前面匹配好的数组，不然达不到加速的目的。对于特定的i, 我们称区间[i, i + z[i]]是i的匹配段，称为Z-box。算法中，我们维护右端点最靠右的匹配端。记作[l, r]。根据定义，s[l, r]是s的前缀。在计算z[i]，我们保证l<=i。初始化时，l=r=0。

def z_function(s):
    n = len(s)
    z = [0] * n 
    l, r = 0, 0 
    for i in range(1, n):
        if i <= r and z[i -l] < r -l + 1:
            z[i] = z[i - l]
        else:
            z[i] = max(0, r - i + 1) # 关键加速步骤
            # 朴素算法
            while i + z[i] < n and s[z[i]] == s[i + z[i]]:
                z[i] += 1
        if i + z[i] - 1 > r:
            l = i
            r = i + z[i] - 1
    return z

复杂度分析: 对于内层的while循环，每次执行都会使得r向后移动至少1位，由于r<n-1所以总共会执行n次。对于外层循环，只有一次线性遍历。总复杂度为O(n)。

应用场景

匹配所有子串

def zMatch(p, t):
    s = p + '$' + t
    Z = z_function(s)
    occurrences = []
    for i in range(len(p) + 1, len(s)):
        if Z[i] >= len(p):
            occurrences.append(i - (len(p) + 1))
    return occurrences

1 2	$ zMatch('abcdabcd', 'abcdefghijklabcdabcd') $ [12]

基因组大小

常说人类基因组大小是3G，指得是3G个碱基。

1 2	$ 3G bp = 10^9 bp $ # 人是二倍体，所以基因组大小为 3G * 2 = 6G

测序数据量

常说的测序数据量是多少Mreads, 多少G碱基数目。 M：用于描述测序reads数，G：用于描述测序获得的碱基总数。

1 2	$ base number = reads number * read length $ 3G bp = 310^9/(1502) = 10M reads # PE150 测序

单位换算

链的方向性

正链和负链基因组中的正链(forward strand)和负链(reverse strand)是有区别的，正链是参考基因组序列，负链是其反向互补序列。
正义链和反义链正义链(sense strand)和反义链(antisense strand)是针对基因转录和翻译过程来定义的；正义链是DNA 双链中携带编码蛋白信息的链，也称为编码链，序列与mRNA序列相同；反义链是与正义链互补的一条链，与mRNA序列反向互补，但是反义链是转录中给mRNA当做模板的链，因此反义链也是模板链。

在基因组中，有的基因的正义链是正链，有的基因的正义链是负链。因此，在测序文库构建中，双链DNA中的一条链对某些基因来说是正义链，对另外一些基因来说是反义链。

总结：图片来源于北大生科院PPT

$ 正义链(sense strand)     == 编码链(coding strand)       == 非模板链(non template strand)
$ 反义链(antisense strand) == 非编码链(non coding strand) == 模板链(template strand)

$ 正链上可以同时存在正义链(sense strand)和反义链(antisense strand), 负链亦然。

链特异性文库

* 普通的RNA-Seq建库方式：第一步，进行RNA到cDNA的反转录，第二步，在反转录以后，普通的RNA-Seq就直接使用random primer进行第2条链合成，随后加adapter，扩增成文库。这样构建出来的RNA-Seq库进行测序以后是分不清这个序列是来自于genome的那条链的，因为被测序的有可能是gene的foward strand，也有可能是reverse strand。

链特异性RNA-Seq建库方式：链特异性建库(以图中间的dUTP方法为例)，首先利用随机引物合成RNA的一条cDNA链(反义链，RNA为正义链)，在合成第二条链的时候用dUTP代替dTTP(正义链)，加adaptor后用UDGase处理，将有U的第二条cDNA降解掉。降解发生之后，双链的文库就只剩下了一条链（反义链）。而这条链的两头是接的不同序列的接头。通过PCR扩增，最终只保留了第一条cDNA（反义链）上机测序。这样最后的insert DNA fragment都是来自于第一条cDNA（反义链），也就是dUTP叫fr-firststrand的原因(测第一条链的fr)。 dUTP测序中的pair-read中的read1(R1)和基因的方向相反，read2(R2)和基因的方向相同。测序是5'->3'进行的。因为read1测的是第一链5'端(反义链)，所以read1和基因的方向相反(反义链)。因为read2测的是第一链的3'端(反义链)，所以，read2和insert序列反向互补(正义链)，和基因的方向相同。因此，dUTP方法是RF。

在reads比对到参考基因组时，那些比对到基因方向（正义链方向）的正链reads就是正义链reads，但是那些比对到基因方向反方向（反义链方向）的正链reads就是反义链reads。那么同样，比对到基因方向的负链reads就是正义链reads，而比对到基因方向反方向（反义链方向）的负链reads就是反义链reads。从而最终将所有正义链reads和反义链reads区分开来。因此在确定基因表达水平时，可以避免基因反义链上的reads匹配的干扰，从而更加准确的检测基因转录表达水平。而且LncRNA的测序也离不开链特异性建库技术。

链特异性文库的优点：1、明确reads链信息的来源，便于基因结构注释；2、更加准确地基因定量

针对链特异性文库的软件设置

RNA-Seq Strand

	forward (transcript)	reverse (rev comp of transcript)
TopHat/Cufflinks `--library-type`	`fr-secondstrand`	`fr-firststrand`
STAR	1st read strand	2nd read strand
Picard CollectRnaSeqMetrics `STRAND_SPECIFICITY`	`FIRST_READ_TRANSCRIPTION_STRAND`	`SECOND_READ_TRANSCRIPTION_STRAND`
htseq-count `-s/--stranded`	`yes`	`reverse`
subread featureCounts `-s`	`1`	`2`
RSEM `--forward-prob`	`1`	`0`
Salmon/Sailfish `--libType`	`SF`/`ISF`	`SR`/`ISR`
HISAT2 `--rna-strandness`	`FR` (`F` for single-end)	`RF` (`R` for single-end)
Library Kit	Illumina ScriptSeq	Illumina TruSeq Stranded Total RNA

基因组坐标系统0-based和1-based

Posted on 2023-12-26 Edited on 2023-12-28

基因组坐标系统

基因组文件各式各样，但是主要呈现的还是序列的相关信息；但是各类注释文件采用的坐标系统却存在差异， 0-based和1-based坐标均被使用。在计算机中，不同计算机语言对字符串索引采用的规则也存在0-based和1-based，因此，选择正确的索引规则会影响着序列处理的准确性。常用的基因组文件格式如下: file format

0-based: 0 1 2 3 4 (UCSC, BED, bedGraph, narrowPeak, BAM)
1-based: 1 2 3 4 (NCBI, Ensembl, GFF, GTF, VCF, SAM, wiggle)
sequence: A T C G

0-based文件第一个碱基前面的间隔为0，第二个间隔为1，依次类推，数字记录的是碱基之间的间隔的顺序。
1-based文件第一个碱基为1，第二个碱基为1，依次类推，数字记录的碱基本身的顺序。

索引

序列$ATCG$中$TC$的索引坐标分别为

$ # 索引坐标
$ 0-based: [1, 3) # 左闭右开
$ 1-based: [2, 3] # 左闭右闭
$ 
$ # 长度计算
$ 0-based: 3 - 1     = 2 # stop - start
$ 1-based: 3 - 2 + 1 = 2 # stop - start + 1

## 计算机语言

不同的计算机语言对字符串的索引不同。

$ # C C++ Python Java 等语言字符串/数组索引采用0-based
$ string = 'ATCG'
$ s1 = string[1:3] # TC 
$ s2 = string[0:2] # AT 
$ # length 
$ len(s1) == 3 - 1 # True
$ len(s2) == 2 - 0 # True

$ # R Matlab Julia COBOL等语言字符串/数组索引采用1-based
$ string <- 'ATCG'
$ # substr(string, start, end)
$ substr(string, 1, 3)   # ATC
$ nchar(s3) == 3 - 1 + 1 # TRUE

0-based/1-based 转化

0-based --> 1-based, start+1, stop保持不变

1-based --> 0-based, start-1, stop保持不变

UCSC

字符串匹配-KMP算法

Posted on 2023-12-13 Edited on 2023-12-21

字符串

子串 substring 字符串中任意个连续的字符组成的子序列称为该字符串的子串(Substring) $S[i..j], i<=j$
子序列 subsequence 字符串$S$的子序列是从$S$中将若干元素提取出来并不改变相对位置形成的序列, $S[p_1], S[p_2], .. S[p_k], 1<=p_1< p_2< p_k <= |S|$。
子串一定是子序列，但是子序列不一定是子串
后缀 suffix
后缀是指从某个位置 i 开始到整个串末尾结束的一个特殊子串, $suffis(S,i) = S[i..|S|-1]$
前缀 prefix
前缀是指从串首开始到某个位置 i 结束的一个特殊子串, $prefix(S, i) = S[0..i], i<=|S|-1$。

字符串匹配就是给定文本字符串(String)是一个长度为$n$的数组$S[1..n]$, 模式(Pattern)是一个长度为$m$且$m<=n$的数组$P[1..m]$，在string中寻找 pattern的过程。

KMP算法

KMP (Knuth-Morris-Pratt) 字符串匹配算法，利用模式串的重复出现的字符实现跳跃，避免主串指针的回溯，仅仅后移模式串。 KMP 算法的核心是PMT(Partial Match Table)部分匹配表数组，PMT的值就是最长公共前后缀的长度，即字符串的前缀组合和后缀组合的交集中最长元素的长度。在模式串中出现前后缀相同时候，在匹配的时候可以实现跳跃，

利用PMT数组减少重复匹配的次数，提升效率。
时间复杂度$O(n+m)$
Next 数组将PTM数组整体右移动一位，原先第一位为-1

第一步，PTM表格的制作


def get_prefix(s: str) -> set:
    if len(s) == 1:
        return set()
    else:
        prefix = set()
        for i in range(1, len(s)):
            prefix.add(s[:i])
        return prefix

def get_suffix(s: str) -> set:
    if len(s) == 1:
        return set()
    else:
        suffix = set()
        for i in range(len(s)):
            suffix.add(s[i:])
        return suffix

def longest_prefix_suffix(suffix: set, prefix: set) -> int:
    inter = suffix & prefix
    if len(inter) == 0:
        return 0
    else:
        return max([len(i) for i in inter])

def getPTM(pat: str) -> list:
    ptm = [0 for i in range(len(pat))]
    for i in range(len(pat)):
        suffix = get_suffix(pat[:i+1])
        prefix = get_prefix(pat[:i+1])
        ptm[i] = longest_prefix_suffix(suffix, prefix)
    return ptm

a = 'abababca'
getPTM(a)
# [0, 0, 1, 2, 3, 4, 0, 1]

# p模式字符串 a b a b a b c a
# 下标        0 1 2 3 4 5 6 7
# PTM         0 0 1 2 3 4 0 1

# PTM[0], p[0]        即"a"的最长公共前后缀，显然为0, 前后缀均不包括字符串本身
# PTM[1], p[0] ~ p[1] 即"ab"的最长公共前后缀，显然为0
# PTM[2], p[0] ~ p[2] 即"aba"的最长公共前后缀，即"a", 长度为1
# PTM[3], p[0] ~ p[3] 即"abab"的最长公共前后缀，即"ab", 长度为2
# PTM[4], p[0] ~ p[4] 即"ababa"的最长公共前后缀，即"aba", 长度为3
# PTM[5], p[0] ~ p[5] 即"ababab"的最长公共前后缀，即"abab", 长度为4
# PTM[6], p[0] ~ p[6] 即"abababc"的最长公共前后缀, 显然为0, 后缀一定含c, 前缀一定含a
# PTM[7], p[0] ~ p[7] 即"abababca"的最长公共前后缀, 即"a", 长度为1,

第二步，Next数组生成


# 根据PTM去求 next 数组：next 数组相当于“PTM” 整体向右移动一位，然后初始值赋为-1。
# Next 数组就可以看成字符串匹配的过程，即，以模式字符串为主串，以模式字符串的前缀为目标字符串，
# 一旦字符串匹配成功，那么Next的值就是匹配成功的字符串的长度
# 这里最核心的点就是j的回溯。
# getNext的时间复杂度是 O(m)
def getNext(p):
    Next = [0 for i in range(len(p))]
    Next[0] = -1 
    i = 0
    j = -1
    while i < len(p)-1:
        # 匹配或者j跳到了最开始
        if j == -1 or p[i] == p[j]:
            i += 1
            j += 1
            Next[i] = j
        # 匹配不成功则在前面的子串中继续搜索
        else:
            j = Next[j]
    return Next

# 若字符串只有一个元素，不执行while，Next[0] = -1
# 如字符串有多个元素,如a='abababca'
# i=0, j=-1; 执行一次while，执行if，i=1,j=0, Next[1] = 0;
# i=1, j=0; 执行while, 如果p[2]=p[0]=a匹配，i=2,j=1, Next[2] = 1; 
#                      如果不匹配，j=Next[0]=-1

a = 'abababca'
getNext(a)
# [-1, 0, 0, 1, 2, 3, 4, 0]

第三步，KMP匹配算法实现

def KMP(pat:str, target:str) -> int:
    i, j = 0, 0
    Next = getNext(pat) # O(m)
    # O(n)
    while i < len(target) and j < len(pat):
        if j == -1 and target[i] == pat[j]:
            i += 1
            j += 1
        else:
            j = Next[j]
    if j == len(pat):
        return i - j
    else:
        return -1 
    
p = "ATG"
t = "ATCGATCGATG"
KMP(t, p)
# 8

字符串匹配-Naive算法

Posted on 2023-12-12 Edited on 2023-12-13

字符串匹配算法

字符串匹配算法在生物信息学中应用广泛，在各种比对算法高频出现。
字符串$S$是将$n$个字符顺序排列而成的序列，$n$为其长度。 - 如果字符串的下标是$1$开始，S的第i个字符是$S[i]$ - 如果字符串的下标从$0$开始，S的第i个字符是$S[i-1]$

字符串匹配就是给定文本字符串(String)是一个长度为$n$的数组$S[1..n]$, 模式(Pattern)是一个长度为$m$且$m<=n$的数组$P[1..m]$，在string中寻找 pattern的过程。

Naive String Match

朴素的字符串匹配算法，也称为暴力匹配算法。

没有预处理
滑动窗口为1
时间复杂度$O((n-m+1)m)$
模式中的字符串比较顺序不限

def naive_match(partten, text):
    occurrences = []
    for i in range(len(text) - len(partten) +1):
        match = True 
        for j in range(len(pattern)):
            if text[i+j] != pattern[j]:
                match = False
                break
        if match:
            occurrences.append(i)
    return occurrences

排序算法

Posted on 2023-12-08 Edited on 2023-12-13

排序算法

排序算法就是将一组数据中的每个元素按照特定的顺序排序。常见的排序算法很多，死记硬背不可取，熟悉原理吾所取。

1. 冒泡算法 bubble sort

循环数组每次选两个数比较，将较小的放置到前面，慢慢地将小的"浮出水面"。冒泡排序的时间复杂度最差为$O(n^2)$，最好为$O(n)$；空间复杂度为$O(1)$。

def bubble_sort(arr: list) -> list:
    # 注意循环的起始和终止位置
    # 如果i=0, 会导致j+1超过数组元素长度
    # 保证 j 和 j+1 能正常索引
    for i in range(1, len(arr)):
        for j in range(0, len(arr)-i):
            # 比较大小
            if arr[j] > arr[j+1]:
                # 浮出水面操作
                arr[j], arr[j+1] = arr[j+1], arr[j]
    return arr

2. 选择排序 selection sort

选择排序, 最直观，直接在数组中找。 - 第一步，在数组中选择最小的元素放置到排序的起始位置， - 第二步，在剩余的元素中继续寻找最小元素，然后放置到已经排序的序列末尾 - 第三步，重复第二步

选择排序的时间复杂度最差$O(n^2)$，最好为$O(n^2)$；空间复杂度为$O(1)$


def select_sort(arr: list) -> list:
    # 首先循环数组，指定第一个为最小索引
    for i in range(len(arr) - 1):
        # 记录最小索引
        minIndex = i
        # 遍历剩下的元素
        for j in range(i, len(arr)):
            if arr[j] < arr[minIndex]:
                minIndex = j
        # 如果i 不是最小值的索引，交换i和minIndex交换
        if i != minIndex:
            arr[i], arr[minIndex] = arr[minIndex], arr[i]
    return arr

3. 插入排序 insertion sort

与打扑克牌的原理类似，它的工作原理为将待排列元素划分为「已排序」和「未排序」两部分，每次从「未排序的」元素中选择一个插入到「已排序的」元素中的正确位置。插入排序的时间复杂度最高为$O(n^2)$，即原始数组倒序排列。时间复杂度最好为$O(n)$，空间复杂度为$O(1)$。

算法稳定性假设在数列中存在a[i]=a[j]，若在排序之前，a[i]在a[j]前面；并且排序之后，a[i]仍然在a[j]前面。则这个排序算法是稳定的！插入排序是稳定的。

def insert_sort(arr: list) -> list:
    # 循环数组
    for i in range(len(arr)):
        preIndex = i - 1 # 前一个索引
        current = arr[i] # 保存当前值
        while preIndex >= 0 and arr[preIndex] > current:
            arr[preIndex+1] = arr[preIndex] # 交换，将前一位的值换到后一位
            preIndex -= 1                   # 索引-1，继续交换；直到换到最前面，preIndex<0停止操作
        arr[preIndex+1] = current           # 当前值的位置
        print(i, arr)
    return arr
# 当i=0时，数组保持原样

# 当i=1时，current = arr[1]
#          preIndex = 0,  如果arr[0] > current, arr[1] = arr[0], preIndex=-1;
# 退出循环,preIndex = -1, arr[0] = curent

# 当i=2时，current = arr[2]
#          preIndex = 1,  如果arr[1] > current, arr[2] = arr[1], preIndex=0; 值交换位置
#          preIndex = 0,  如果arr[0] > current, arr[1] = arr[0], preIndex=-1;
# 退出循环,preIndex = -1, arr[0] = current

# 当i=3时，current = arr[3]
#          preIndex = 2,  如果arr[2] > current, arr[3] = arr[2], preIndex=1;
#          preIndex = 1,  如果arr[1] > current, arr[2] = arr[1], preIndex=0;
#          preIndex = 0,  如果arr[1] > current, arr[1] = arr[0], preIndex=-1;
# 退出循环,preIndex = -1, arr[0] = current     

# 注意：current 是传递的中间值；传递是在相邻元素之间交替进行的。

4. 希尔排序 shell sort

希尔排序是对插入排序的一种优化，插入排序的原理可知，每次只能移动一个位置。 shell sort采用步长分割后再进行排序，本质是分组插入排序。具体实现是，首先数组长度除以步长(gap)，后续步长再除以gap的缩小值，直到以gap=1为步长排序。

def shell_sort(arr: list) -> list:
    import math 
    # 指定步长的长度
    gap = 1
    while (gap < len(arr)/3):
        gap = gap*3 + 1
        print('gap=',gap)
    # 步长长度随着数组的长度而改变，
    # 但是进入排序之前会是一个定值
    while gap > 0:
        #使用直接的插入排序
        for i in range(gap, len(arr)):
            temp = arr[i]
            j = i - gap 
            while j >= 0 and arr[j] > temp:
                arr[j+gap] = arr[j]
                j -= gap 
            arr[j+gap] = temp
        # 减小gap
        gap = math.floor(gap/3)
    return arr

测序碱基质量值

Posted on 2023-12-07 Edited on 2023-12-13

测序碱基质量值

对于高通量测序来说，碱基质量值是衡量测序数据的重要指标。碱基质量值常用Phred Quality Score 表示，简写Q。 p值表示碱基被测错误的概率, p越小，Q越大。

\[ Q=-10log_{10}p \]

Q&p转换

\[ p = 10^{Q \over -10} \]

python Q-p转换实现

def error_prob(quality):
    qval = quality 
    return 10**(qval/-10.0)

print("{0:^5} {1:^10}".format("Phred", "Prob of"))
print("{0:^5} {1:^10}".format("score", "error"))
for phred in range(0,42):
    print("{0:^5} {1:03.5f}".format(phred, error_prob(phred)))

Phred  Prob of  
score   error   
  0   1.00000
  1   0.79433
  2   0.63096
  3   0.50119
  4   0.39811
  5   0.31623
  6   0.25119
  7   0.19953
  8   0.15849
  9   0.12589
 10   0.10000
 11   0.07943
 12   0.06310
 13   0.05012
 14   0.03981
 15   0.03162
 16   0.02512
 17   0.01995
 18   0.01585
 19   0.01259
 20   0.01000
 21   0.00794
 22   0.00631
 23   0.00501
 24   0.00398
 25   0.00316
 26   0.00251
 27   0.00200
 28   0.00158
 29   0.00126
 30   0.00100
 31   0.00079
 32   0.00063
 33   0.00050
 34   0.00040
 35   0.00032
 36   0.00025
 37   0.00020
 38   0.00016
 39   0.00013
 40   0.00010
 41   0.00008

Ascii Codes

测序文件中使用的是Ascii Codes字符，每个字符代表一个phred值。 Ascii Codes字符在计算机中和整型数字一一对应。

print ("{0:^8}  {1:^8}".format("Character", "ASCII #"))
for i in range(33,90):
    print("{0:^8}  {1:^8}".format(chr(i),i))
Character  ASCII #
   !         33   
   "         34   
   #         35   
   $         36   
   %         37   
   &         38   
   '         39   
   (         40   
   )         41   
   *         42   
   +         43   
   ,         44   
   -         45   
   .         46   
   /         47   
   0         48   
   1         49   
   2         50   
   3         51   
   4         52   
   5         53   
   6         54   
   7         55   
   8         56   
   9         57   
   :         58   
   ;         59   
   <         60   
   =         61   
   >         62   
   ?         63   
   @         64   
   A         65   
   B         66   
   C         67   
   D         68   
   E         69   
   F         70   
   G         71   
   H         72   
   I         73   
   J         74   
   K         75   
   L         76   
   M         77   
   N         78   
   O         79   
   P         80   
   Q         81   
   R         82   
   S         83   
   T         84   
   U         85   
   V         86   
   W         87   
   X         88   
   Y         89

Phred 编码

Phred+33
phred+33编码，就是使用phred 值 + 33,求和之后对应的ascii字符就是该碱基的字符编码。
1
2
$ Phred Sum Ascii
$ 30 63(30+33) ?
Phred+64 phred+64定义类似，就是使用phred 值 + 64,求和之后对应的ascii字符就是该碱基的字符编码。比较老的测序数据采用该编码方式。
1
2
$ Phred Sum Ascii
$ 30 94(30+64) ^
Phred、Error、Phred+33和Phred+64

print ("{0:^5}  {1:^8}  {2:^8}  {3:^8}".format("Phred",  "Prob of", "Phred+33", "Phred+64"))
print ("{0:^5}  {1:^8}  {2:^8}  {3:^8}".format("score",  "Error", "Ascii", "Ascii"))
for phred in range(0,42):
    print ("{0:^5}  {1:03.5f}  {2:^8}  {3:^8}".format(phred, error_prob(phred), chr(phred+33), chr(phred+64)))

Phred  Prob of   Phred+33  Phred+64
score   Error     Ascii     Ascii
  0    1.00000     !         @
  1    0.79433     "         A
  2    0.63096     #         B
  3    0.50119     $         C
  4    0.39811     %         D
  5    0.31623     &         E
  6    0.25119     '         F
  7    0.19953     (         G
  8    0.15849     )         H
  9    0.12589     *         I
 10    0.10000     +         J
 11    0.07943     ,         K
 12    0.06310     -         L
 13    0.05012     .         M
 14    0.03981     /         N
 15    0.03162     0         O
 16    0.02512     1         P
 17    0.01995     2         Q
 18    0.01585     3         R
 19    0.01259     4         S
 20    0.01000     5         T
 21    0.00794     6         U
 22    0.00631     7         V
 23    0.00501     8         W
 24    0.00398     9         X
 25    0.00316     :         Y
 26    0.00251     ;         Z
 27    0.00200     <         [
 28    0.00158     =         \
 29    0.00126     >         ]
 30    0.00100     ?         ^
 31    0.00079     @         _
 32    0.00063     A         `
 33    0.00050     B         a
 34    0.00040     C         b
 35    0.00032     D         c
 36    0.00025     E         d
 37    0.00020     F         e
 38    0.00016     G         f
 39    0.00013     H         g
 40    0.00010     I         h
 41    0.00008     J         i

为什么是+33

因为在计算机中，0-32整型对应的字符有很多是空白和非打印的字符，肉眼不可见，所以才在33位开始编码。

1
2
3

print ("{0:^8}  {1:^8}".format("Character", "ASCII #"))
for i in range(0,40):
    print("{0:^8}  {1:^8}".format(chr(i),i))