生信编程系列（1-2）
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
1.生信编程之思维讲解
有很多的生信软件可以帮助我们处理分析生物学数据，但在没有这些软件的时候，我们可以自己编写一些程序来实现相同的目的，这些程序可能是一个很长的脚本，也可能只是一条简单的shell命令，自己学会编写一些脚本可以更好地理解数据的处理分析过程。
对FASTQ文件进行操作
fastqc可以统计出fastq文件中的序列总数
如果要自己通过编程来统计总的序列总数，我们首先要了解fastq格式的意义。
fastq文件是一种存储来生物序列（通常是核酸序列）以及相应质量评价的文本格式 fastq每四行代表一条reads的记录 第一行通常以'@’开头，后面紧跟着ID以及其他信息，包括测序设备以及样本编号等。 第二行记录的是read序列 第三行以'+’开头，后面是reads的名称（一般与@后面的内容相同） 第四行代表reads的质量。起初Sanger中心使用Phred quality score来衡量read中每个碱基的质量，即-10lgP，其中P代表的是该碱基被测序错误的概率，如果该碱基测序出错的概率为0.001，Q应该为30，那么30+33=63，63对应的ASCii码为“？”
了解了fastq格式之后，我们要想统计一个fastq文件中记录的reads总数就可以通过一行简单的命令实现啦
wc -l longreads.fq
==注意== 上面这了shell命令统计的是fastq文件中全部的行数，但我们已经知道fastq中每4行为一个记录，那么就用该命令计算出来的值除以4，就是我们需要的值，假设该命令计算的结果为24000，那么总共的reads数就是6000.
脚本对于生信人来说是非常重要的，尤其在处理很大的数据的时候
截掉5'端和3'端的几个碱基
由于二代测序技术的缺陷，测序结果的3'端和5’端的碱基质量相对较差，为了避免这些低质量的碱基下游分析产生不必要的影响，我们必须将3'端和5’端一些低质量的碱基去掉。虽然很多软件都可以帮我们实现这个目的，但是我们接下来将通过自己编写命令来实现。
我们已经知道fastq文件的第二行为reads序列，第四行为对应的质量值，我们不能只把碱基去掉，而不考虑质量值。
perl -alne '{print}' longreads.fq | head
==解释== 该命令的作用是将文件打印出来，然后用head查看前10行
perl -alne '{print if $.%4==1}' SRR2188201.1_2.fastq | head
==解释== perl中 $.代表的是行，每4行取一行，因为fastq文件的总行数为4的倍数，那么以四行为一个单位，四行中的第一行所在位置的行数除以4，余数为1，所以用$.%4==1返回的就是reads记录的第一行。
perl -alne '{print length}' SRR2188201.1_2.fastq | head
==解释== 返回每一行的长度
perl -alne '{print if length($_) > 100}' SRR2188201.1_2.fastq | head
==解释== 因为含有reads序列和reads碱基质量这两行的长度是一致的，而且相比于另外两行来说长度也不一样，所以可以根据长度筛选出第二行和第四行
perl -alne '{print length if $.%4==2}' SRR2188201.1_2.fastq | head
==解释== 每四行中只返回第二行的长度
对第二列的reads序列和第四列的碱基质量进行修剪，3'端和5’端分别去掉4个碱基
perl -alne '{print substr($_,9) if $.%4==2|$.%4==0;print $_ if $.%4==1|$.%4==3}' SRR2188201.1_2.fastq | head
==解释== 我们只需要对FASTQ文件中的每四行中的第二行和第四行进行处理，另外两行不需要处理的，就原样输出
序列长度分布统计
perl -alne '{print length if $.%4==2}' SRR2188201.1_2.fastq | perl -alne '{$h{$_}++}END{print "$_\t$h{$_}" foreach sort{$a <=> $b} keys %h}'
==解释== 首先返回四行记录中的第二行，也就是reads序列，使用哈希进行统计，如果长度相同，数量就加1。
perl -alne '{print length if $.%4==2}' SRR2188201.1_2.fastq | perl -alne '{$tmp=int($_/10);$h{$tmp}++}END{print "$_\t$h{$_}" foreach sort{$a <=> $b} keys %h}'
规定以区间长度为10来进行统计条数
FASTQ 转换成 FASTA
首先我们要了解FASTA格式和FASTQ格式的差别，FASTA格式是一种基于文本用于表示核酸序列或者多肽序列的格式。其中核酸或者氨基酸均以单个字母来表示。FASTA格式没两行记录一条序列，第一行以'>'符号开头，后面紧跟着序列名称，第二行记录的是序列。FASTQ格式相比于FASTA第三行的名称和第四行的碱基质量，其次是第一行名称的开头符号不一致。
perl -alne '{print $_ if $.%4==1|$.%4==2}' SRR2188201.1_2.fastq | perl -alne '{s/@/>/g;print}' > SRR2188201.1_2.fasta
==解释== 首先，通过用行号处以4判断目前是第几行，将第三第四行去掉，只留下第一行序列名称和第二行碱基序列。通过观察可以发现，第一行序列名称的开头均为@，而且在其他位置没有@出现，所以就可以使用替换的方式，加将全部的@替换为>
统计碱基个数及GC%
perl -alne '{print $_ if $.%4==2}' SRR2188201.1_2.fastq | perl -alne '{@CG=($_=~m/(G|C)/g);$CG=@CG; $sum_CG+=$CG; $sum_all+=length($_)}END{print $sum_CG/$sum_all}'
##对FASTA文件的操作
FASTA文件格式是一种基于基于文本用于表示核酸序列或者多肽序列的格式，FASTA文件中每两行数据记录一条reads，其中第一行一般以'>'开头，后面紧跟着序列名称，第二行是核酸序列或者多肽序列。
取出反向互补的序列
perl -alne '{print $_ if $_=~/(.)(.)(.)(.)(.)\5\4\3\2\1/;print $&}' SRR2188201.1_2.fasta
==解释== 是痛perl正则表达可以匹配回文序列
取序列的互补序列
perl -alne '{print $_ if $.%2==1;$com=$_ if $.%2==0;$com=~tr/tcgaTCGA/agctAGCT/;print $com}' SRR2188201.1_2.fasta | head
==解释== 通过判断行当前行是否为包含序列的行，将序列名称先输出，然后对序列进行碱基替换，根据A-T， C-G配对的原则进行替换
###取反向序列
perl -alne '{print $_ if $.%2==1; print $com = reverse($_) if $.%2==0}' SRR2188201.1_2.fasta | head
==解释== 使用reverse函数就可以得到反向序列
获得反向互补序列
perl -alne '{print $_ if $.%2==1;$com=$_;$com=~tr/tcgaTCGA/agctAGCT/;print reverse($com) if $.%2==0}' SRR2188201.1_2.fasta | head
DNA to RNA
perl -alne '{print $_ if $.%2==1; $com=$_; $com=~tr/tT/uU/;print $com if $.%2==0}' SRR2188201.1_2.fasta | head
==解释== 要将DNA转变为RNA，就需要将DNA中的碱基T替换为碱基U
模拟DNA转录为RNA
perl -alne '{print $_ if $.%2==1;$com=$_;$com=~tr/atcgATCG/uagcUAGC/;print reverse($com) if $.%2==0}' SRR2188201.1_2.fasta | head
==解释== DNA转录为RNA，主要基于碱基配对，而且RNA和模版DNA互为反向互补序列
大小写字母形式输出
#将核酸序列以小写输出
perl -alne '{print $_ if $.%2==1; print lc($_) if $.%2==0}' SRR2188201.1_2.fasta | head#将核酸序列以大写输出
perl -alne '{print $_ if $.%2==1; print uc($_) if $.%2==0}' SRR2188201.1_2.fasta | head
==解释== perl函数uc() [uppercase]，lc()[lowercase]能够进行大小写转换
###每行执行长度输出序列
使用python来进行处理
#!/usr/bin/python
import argparse
import re
parser = argparse.ArgumentParser()
parser.add_argument('infile',help='input genome file of fasta',type=str)
parser.add_argument('width',help='the number of bases in a output line',type=int)
args = parser.parse_args()
infile = args.infile
width = args.width
seq_base={}
def get_chr(buf):
(buf,tmp) = buf.split('>',1)
if '>' in tmp:
get_chr(tmp)
return buf,tmp
with open(infile,'r') as f:
tmp = f.readline()
#print(tmp)
chr_id = re.split(r'\s',tmp)[0][1:]
#print(chr_id)
seq_base[chr_id] = []
while 1:
buffer = f.read(1024*1024)
if not buffer:
break
if '>' in buffer:
(buffer,tmp) = get_chr(buffer)
buffer = buffer.upper()
seq_base[chr_id].append(buffer.replace('\n',''))
(tmp,buffer) = tmp.split('\n',1)
chr_id = re.split(r'\s',tmp)[0]
seq_base[chr_id]=[]
buffer = buffer.upper()
seq_base[chr_id].append(buffer.replace('\n',''))
def seq_width(seq,width):
start=0;end=start+width
while end < len(seq):
print(seq[start:end])
start+=width;end+=width
while end >= len(seq):
print(seq[start:])
for seqname,seq in seq_base.items():
if seqname == 'chr1':
seq = ''.join(seq)
seq_width(seq,width)
==解释== 首先创建一个字典来存储序列信息，结构如下
{'chr1':'ATGCTAGCTAG...',
'chr2':'ATGCTAGCTAG...'}
然后创建一个函数按照规定的碱基数量输出
def seq_width(seq,width):
start=0;end=start+width
while end < len(seq):
print(seq[start:end])
start+=width;end+=width
while end >= len(seq):
print(seq[start:])
注意 以下步骤都是基于上述脚本创建的序列字典进行的
###按照序列长度/名字排序
按照序列名称进行排序
for i in sorted(seq_base):
print(i+'\n'+seq_base[i])
seq_base是一个存储了基因组序列名称和序列信息的字典，键为序列名称，值为序列信息，sorted(dict)可以直接根据键对字典进行排序，对于 seq_base来说就是根据序列名称进行排序
按照序列长度进行排序
len_base = {}
for seqname,seq in seq_base.items():
len_base[seqname]=len(''.join(seq))
for name,seq in sorted(len_base.items(),key = lambda kv:(kv[1],kv[0])):
print(name,seq)
在之前包含序列的那个字典的基础上，重新创建一个字典，字典的键为序列名称，而字典的值为序列的长度。通过创建一个匿名函数来实现根据值对字典进行排序，在len_base这个字典中，值就是序列的长度。sorted(chrs,key=lambda)可以根据自定义的顺序进行排序。
提取指定ID的序列
for seqname,seq in seq_base.items():
if seqname == "chr1":
print('>'+seqname+'\n'+''.join(seq))
==解析== seq_base是一个键为序列名称，值为序列的字典，通过字典查询可以获得目的ID序列
随机抽取序列
random_chr = random.sample(seq_base.keys(),1)[0]
print('>'+random_chr+'\n'+''.join(seq_base[random_chr]))
==解析== random.sample()可以从列表中随机取出指定的元素，再以列表的形式返回，所以要用[]取出目的元素
高级难度
根据坐标取序列
首先使用bedtools取出序列，作为参照
bedtools getfasta -fi hg38.fa -bed fasta.bed -fo tmp.fa
==参数介绍== -fi输入的fasta格式文件,-bed包含所需要截取序列的染色体号和序列起始和终止位置,-fo输出文件fasta.bed文件：
chr1 43045448 43046048
chr2 43045461 43046061
chr3 43056769 43057369
chr4 43056816 43057416
使用python脚本进行提取
#!/usr/bin/python
import pysam
hg38 = pysam.FastaFile('hg38.fa')
with open('fasta.bed','r') as f:
for i in f.readlines():
if i.strip() != '':
(chrom,start,end) = i.strip().split('\t')
print(chrom+':'+'\t'+start+'-'+end+'\n'+hg38.fetch(chrom)[int(start):int(end)])
最终得到的结果于bedtools取出的一样，最主要的是用到pysam这个模块
2.生信编程之hg38基因组序列探究
统计hg38基因组中的GC含量
统计N和GC含量==> 分别统计ATGCN的数量，然后使用（N数量/总数，G数量+C数量/总数）
首先要构建数据结构，也就是说用这种数据结构来统计碱基数量 python中的字典
{
'chr1':{
'A':n,
'T':n,
'C':n,
'G':n,
'N':n
},
'chr2':{
'A':n,
'T':n,
'C':n,
'G':n,
'N':n
},
}
初级脚本
hg38 = open('test.fa','r')
sum_atgc = {}
for line in hg38:
if line.startswith('>'):
chr_id = line[1:]
sum_atgc[chr_id] = {}
sum_atgc[chr_id]['A'] = 0
sum_atgc[chr_id]['T'] = 0
sum_atgc[chr_id]['G'] = 0
sum_atgc[chr_id]['C'] = 0
sum_atgc[chr_id]['N'] = 0
else:
line = line.upper()
sum_atgc[chr_id]['A'] += line.count('A')
sum_atgc[chr_id]['T'] += line.count('T')
sum_atgc[chr_id]['G'] += line.count('G')
sum_atgc[chr_id]['C'] += line.count('C')
sum_atgc[chr_id]['N'] += line.count('N')
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
print("A : %s" % (atgc_count['A']))
print( "T : %s" % (atgc_count['T']))
print( "G : %s" % (atgc_count['G']))
print( "C : %s" % (atgc_count['C']))
print( "N : %s" % (atgc_count['N']))
print( "total : %s" % (sum(atgc_count.values())))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
虽然逻辑上符合，而且也能够计算出目的的结果，但是重复代码太多，需要进行优化 ==注意== 由于基因组中可能同时存在大小写不同的碱基，所以在进行统计地时候需要统一大小写line.upper()
减少代码重复问题
#!/usr/bin/python
import sys
args = sys.argv
hg38 = open(args[1],'r')
bases = ['A','T','G','C','N']
sum_atgc = {}
for line in hg38:
if line.startswith('>'):
chr_id = line[1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
else:
line = line.upper()
for base in bases:
sum_atgc[chr_id][base] += line.count(base)
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
for base in bases:
print("%s : %s" % (base, atgc_count[base]))
print( "total : %s" % (total))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
==注意== 将需要进行处理的碱基放到一个列表中，然后对列表进行循环，这样就能避免重复写相同的代码
从1开始写python
0.5 开始课程 1 基础语法 2 编程技巧（打开文件，参数处理等） 3 效率（程序的运行时间）
提高效率
#!/usr/bin/python
import sys
import time
import re
start = time.process_time()
args = sys.argv
sum_atgc = {}
bases = ['A','T','G','C','N']
def get_chr(buf):
(buf, tmp) = buf.split('>',1)
if '>' in tmp:
get_chr(tmp)
return buf,tmp
with open(args[1],'r') as f:
tmp = f.readline()
chr_id = re.split(r'\s',tmp)[0][1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
while 1:
buffer = f.read(1024 * 1024)
if not buffer:
break
if '>' in buffer:
#print(len(buffer.split('>')))
(buffer,tmp) = get_chr(buffer)
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
(tmp, buffer) = tmp.split('\n', 1)
chr_id = re.split(r'\s', tmp)[0][1:]
sum_atgc[chr_id] = {}
for base in bases:
sum_atgc[chr_id][base] = 0
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
else:
buffer = buffer.upper()
for base in bases:
sum_atgc[chr_id][base] += buffer.count(base)
for chr_id,atgc_count in sum_atgc.items():
print(chr_id)
GC = int(atgc_count['G']) + int(atgc_count['C'])
total = sum(atgc_count.values())
for base in bases:
print("%s : %s" % (base, atgc_count[base]))
print( "total : %s" % (total))
print("CG_per : %s" % (GC/total))
print("N_per : %s" % (int(atgc_count['N'])/total))
end = time.process_time()
print("time : %s" % (int(end-start)))
==解析1== 使用f.read()读取文件比较慢，通过使用buffer读取可以提升速度，每次读取的字节数为（1024*1024）。对于fasta格式的基因组文件来说，需要区分'>'开头的序列名和序列。首先要读第一行的序列名称tmp=f.readline()，后面的内容就可以按照buffer的要求来读取。一个buffer读入的文件有三种形式：
ATCGATGACTGA\n >chr1 \nGCTGATGAC
包含序列，中间混入一个序列名称 ATCGATGCA\n> chrx \nAGCTGATCGATCGATCGATGATCGATGCATGCTA\n >chry \nATGCTAGCTAGTC
对于某些比较短的染色体来说，可能存在序列中混入了两个序列名称的情况 ATCGATGCGGCTAGCTAGCTGATCGAGGCCGATGACG
只有序列信息，没有序列名称
==解析2== 使用while循环读入buffer时，一定要判断是否有buffer，否则文件读完之后，也不会退出循环
==解析3== 对buffer中含有两个甚至多个以'>'开头的序列名的情况进行处理：最好的解决方式是递归
def get_chr(buf):
(buf, tmp) = buf.split('>',1)
if '>' in tmp:
get_chr(tmp)
return buf,tmp