手记

perl 命令行实战2 - fastq文件的相关操作

慕码人8056858
2018-12-29 23:03:40浏览 3575

1. fastq文件的介绍

详细介绍：fastq格式

fastq是一种除了fasta文件之外，做生信另外一个最为常见的文件类型。

fastq文件主要是用来存储测序的序列信息的。它以4行为一个单位，是严格按照行限定的格式，不像之前的fasta文件那样比较灵活的形式。
虽然网上有很多关于fastq文件的介绍，这里我还是想啰嗦几句，因为这样到后面写perl程序的时候可以前后对应着看会稍微方便一点

示例

@SRR2177462.1 FCC600JACXX:2:1101:1489:2045/1NGGCAAAAGGAAGCACATATTCGCATATAGAACCAGGATTTATAAGGTACAACAANTAGACTTATCCTCCACTCTCATGTTCATGAATCA
+#1=ABDD?FH?HFG>DGHBHIFFHGCGGCHGCGHGII)?FG@DBAEH9???FGAB#-5@(@=EEHECAAH@EE;BCEEFA@ADDCCA;AC@SRR2177462.2 FCC600JACXX:2:1101:1661:2085/1NAATGAAATTAAAGATAGCTGATCTATATTTCTCAAGTGACTAAGTATTAATATTATGCGTACTCTGTATTTCTCTAGTTGGTGGTTTAG
+#4=DDFFFHHFHFIGIHIIIJEGHIGIAHIJFIEHHHIIIIIEIHIIIJCFIHJIHIIAHIFFGIIJHFHIIJIJF@FEHGIDHHA?B@C

现在我们从头开始，每4行为一组，在每一组之内的：

@SRR2177462.1 FCC600JACXX:2:1101:1489:2045/1

  • 第一行：以@字符开头，表示这条序列的测序信息，来自于哪台测序仪、第几个run、第几个lane、第几个tail、以及在这个tail的X轴与Y轴的位置。

NGGCAAAAGGAAGCACATATTCGCATATAGAACCAGGATTTATAAGGTACAACAANTAGACTTATCCTCCACTCTCATGTTCATGAATCA

  • 第二行：测序得到的序列。也可以说是read

+

  • 第三行：以+字符开头，后面跟着ID序列标识符，如果+后面有东西那么必须与第一行一致，一般都是没有的，为了节省内存。

#1=ABDD?FH?HFG>DGHBHIFFHGCGGCHGCGHGII)?FG@DBAEH9???FGAB#-5@(@=EEHECAAH@EE;BCEEFA@ADDCCA;AC

  • 第四行：是第二行序列每个碱基的质量编码值，包含与第二行字符数量相同的符号。经过转换之后可以得到这个碱基的质量数值，以及这个碱基的错误率。

前期准备

因为一般的fastq文件都是经过压缩的，如果你想拿fastq的实际文件来做这次的测试，那样可能比较耗时间，特别是在gzip或者pigz的解压时间上面。所以这里如果你想快速进行实战演练的话，我建议先把本文末尾的6. 测试数据这一步做了，得到123.fastq.gz进行测试会快一些。

注意!!

由于我演示是使用的windows系统，在文件的一行结束的位置除了一个\n换行符之外，还有\r回车符这样的字符存在，而使用perl中的chomp方法不能除去\r回车符，所以下面代码中，在Mac或者Linux中可以直接写为chomp我换成了s/\r?\n//。在后面生成123.fastq.gz文件的时候粘贴复制到文本里面去的时候也可能带上回车符\r，这里使用windows进行测试要注意了。

2. 获取fastq文件的信息

2.1 序列条数

思路 ： 利用fastq的四个为一个单位的特征

  • 方法1
    比较丑陋

zcat 123.fastq.gz | perl -n -e '
  $n++;
  # 因为四个为一组，没必要将所有的行都加起来再除以4
  # 消耗掉其他三行
  <>;<>;<>;
  END{print $n}
'

  • 方法1.0
    也可以这么写

zcat 123.fastq.gz | perl -n -e '
  if($. % 4 == 1){
    $n++;
  }
  END{print $n}
'

或者使用pigz，pigz要比gzip解压快一些，但是git for windows上没有这个东西，另外使用zcat进行查看文件我发现比直接解压文件要快，所以后面测试还是使用zcat。不过也展示一下pigz的使用方式：

# -d : 解压模式# -c : 将解压的得到的数据输出到标准输出中pigz -d -c 123.fastq.gz | perl -n -e '
  $n++;
  <>;<>;<>;
  END{print $n}
'

  • 方法2
    使用linux命令

# wc -l 计算行数line=`zcat 123.fastq.gz | wc -l`let seq_num=$line/4echo $seq_num

输出为

20

2.2 质量值的类型

2.2.0 有关质量值的编码方式

现在一般都是采用的Phred+33类型的，那为什么还要说这一小节呢？

通过介绍这些质量值类型可以帮助理解fastq中质量值与序列的关联性以及在这种利用文本传达信息的过程中，人们所展现出来的智慧！

2.2.0.1 测序的错误来源
+==============+       +==============+| .   .   .   .|       | *   *   *   *||   .   .   .  |  ==>  |   *   *   *  || .   .   .   .|       | *   *   *   *||   .   .   .  |       |   *   *   *  |+==============+       +==============+
 测序芯片中测序点        激发之后产生荧光

在测序的时候，是通过对应位置发出来的对应颜色的荧光的强度来得到这个位置测定出。


chromatogram样图.jpg

那测序的错误是从何而来呢？

在荧光染料测序中，每次发生碱基合成时会释放出荧光信号，该信号被CCD图像传感器捕获。记录下荧光信号的峰值，生成一个实时的轨迹数据（chromatogram）（如上图所示）。因为不同的碱基用不用的颜色标记，检测这些峰值即可判断出对应的碱基。但由于这些信号的波峰、密度、形状和位置等是不连续或模糊的，有时很难根据波峰判断出正确的碱基。通过对每一个点对不同颜色的曲线参与的程度进行计算，计算许多与波峰大小和分辨率相关的参数，根据这些参数，从一个巨大的查询表中找出碱基质量得分。这个查询表是根据对已知序列的测序数据分析得到的（应该是分析得到波峰参数与碱基错误率的关系，再通过公式把错误率转换成质量得分，得到波峰参数与质量得分的直接对应表）。


波长1波长2波长3波长4
A0.650.140.010.02
B0.811.030.020.03
C0.020.020.950.04
D0.030.031.101.56

表2.2.0.1 例如一个测序位点不同碱基所带的荧光基团在这个点发出的荧光的贡献

那么也就是说，一个点ATGC四个碱基上经过激发之后得到不同的波长，如果这个位点不同颜色占据的量越多，那么出错的可能性就越大，如果某种颜色占据的越多那么越有可能是这个碱基。那么如何描述这种碱基的纯粹程度呢？

2.2.0.2 错误的评判标准

Quality Score —— 评判一个碱基判读可靠性指标

Q = -10 lg P

  • Q是质量值

  • P是错误率

准确度P值(错误率)Q值（质量值）
90%10%10
99%1%20
99.9%0.1%30
99.99%0.01%40
99.999%0.001%50

平常常常说Q30，就是说碱基准确度达到99.9%的程度，就是说出错的可能性为0.1%。这里别人推荐了一种记忆方法，10是一个9，20两个9，40是4个9。

这里有关的测序步骤以及质量值可以参看陈老师的视频

陈巍学基因

那么就需要对每次每个测得的碱基的质量进行记录，怎么来记录呢？可是质量值从个位数到十位数会发生什么呢？

2.2.0.3 质量值的记录方式

使用数值的序列与质量值得对应关系

本来对应关系是这样的

N 2
A 19
A 28
T 35
G 35
A 37
A 37
A 37
T 39
T 39
A 37
A 39
A 37
G 40
A 38
. .
. .
. .

写到一行上（其实没有这种数值类型的表示方式）

序列 NAATGAAATTAAAGATAGCTGATCTATATTTCTCAAGTGACTAAGTATTAATATTATGCGTACTCTGTATTTCTCTAGTTGGTGGTTTAG
质量 
21928353537373739393739374038403940404041363839403840323940413740363939394040404040364039404040413437403941403940403239403737384040413937394040414041373137363938403539393230333134

就单单只看两个碱基，碱基与质量值之间的对应关系有多种

N => 2A => 1或者
N => 21A => 9或者
N => 2A => 19或者
N => 21A => 92

如果这么写的话，第一个碱基质量值为个位数，其他的为十位数，这样记录之后，后续程序怎么来分析它呢？一个碱基对应的质量值到底是个位数还是十位数呢？这个样子质量值与碱基之间怎么对应啊？这是个问题......

那如果用空格将数值间隔开呢？

序列 NAATGAAATTAAAGATAGCTGATCTATATTTCTCAAGTGACTAAGTATTAATATTATGCGTACTCTGTATTTCTCTAGTTGGTGGTTTAG
质量 2 19 28 35 35 37 37 37 39 39 37 39 37 40 38 40 39 40 40 40 41 36 38 39 40 38 40 32 39 40 41 37 40 36 39 39 39 40 40 40 40 40 36 40 39 40 40 40 41 34 37 40 39 41 40 39 40 40 32 39 40 37 37 38 40 40 41 39 37 39 40 40 41 40 41 37 31 37 36 39 38 40 35 39 39 32 30 33 31 34

有没有感觉到瞬间多了好多字符一样的，这样带来的后果就是占用的内存急剧增加。

于是乎人们就想到了一个好办法


ASCII码表.jpg

使用ASCII码的序列与质量值得对应关系

序列 NAATGAAATTAAAGATAGCTGATCTATATTTCTCAAGTGACTAAGTATTAATATTATGCGTACTCTGTATTTCTCTAGTTGGTGGTTTAG
质量 #4=DDFFFHHFHFIGIHIIIJEGHIGIAHIJFIEHHHIIIIIEIHIIIJCFIHJIHIIAHIFFGIIJHFHIIJIJF@FEHGIDHHA?B@C

可以看到一个碱基对应一个表示质量值的字符，这样既减少了文件所占据的内存，也让碱基质量值与碱基的对应关系更加清晰。

那么将使用ASCII码进行记录的质量值与直接使用数字进行记录的质量进行比较，使用数字表示法：

  • 文件占据内存增加：在fastq文件中，占内存的大头是序列以及质量值的文字，如果质量值用两个数字表示那内存占据会增加很多

  • 碱基与质量值对应的关系：如果是以数字的形式对应碱基的话，如果一直是1个碱基对应2个数字那都好说，可是保不齐某些碱基测序质量极差，质量值为个位数，那么既有1对2、也有1对1，怎么才能确定碱基与质量值之间的对应关系呢？

既然采用ASCII码来记录碱基质量，那么为什么还要分什么Phred+33、Phred+64之类的？

2.2.0.4 质量值的编码方式

现有的常见的测序质量值编码方式

  • Phred+33

Quality encoding: !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI
                  |         |         |         |         |
   Quality score: 0........10........20........30........40

  • Phred+64

Quality encoding: @ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefgh
                  |         |         |         |         |
   Quality score: 0........10........20........30........40

  • Solexa+64

Quality encoding: ;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefgh
                  |    |         |         |         |         |
   Quality score:-5....0........10........20........30........40

话说Phred和solexa是什么意思呢？百度Phred之后发现有一个很古老的网站，网页比较简洁，里面写着

Phred is a base-calling program for DNA sequence traces. Phred reads DNA sequence chromatogram files and analyzes the peaks to call bases, assigning quality scores ("Phred scores") to each base call.

就是说Phred这个原来是一个工具，用来对测序得到的chromatogram文件分析（可以见上面的有四条曲线的图）。这个软件通过对四种不同的荧光的强度来分析得到该碱基位点的质量值。

当时是1998年发现核酸研究上的一篇文章《 Base-calling of automated sequencer traces using phred. I. Accuracy assessment》然后同年又发了文章《Base-calling of automated sequencer traces using phred. II. Error probabilities》。

最初开发时，Phred在所检查的数据集中产生的错误明显少于其他方法，平均减少40-50％的错误。Phred质量评分已被广泛接受用于表征DNA序列的质量，并且可用于比较不同测序方法的功效 。Phred在人类基因组计划中发挥了重要作用，自动脚本处理大量序列数据。它当时是学术和商业DNA测序实验室最广泛使用的碱基调用软件程序，因为它具有很高的碱基调用准确性。 ——wikipedia

而solexa是一种基于边合成边测序的技术。通过单分子阵列实现小型阵列上的桥式PCR反应。新的阻断技术能够每次只合成一个碱基，经过激光激发之后得到激发光，对激发光进行分析得到碱基信息。

早些年，人们在谈到新一代测序仪时，经常会提起Solexa，而不是Illumina。这是一家低调的公司，规模也不大，但是测序技术却非常新颖。它开发出的测序仪，在通量上领先于其他竞争产品。收购Solexa，也成为Illumina的转折点，从此踏上高速发展的道路。 ——生物通

不同编码类型.png

现在大部分都是用的Illumina 1.8+的Phred+33方式，也就是将质量值加上33，得到一个数值，再将这个数值按照上面的ASCII码表对应到特定的字符上去。

比如

Q值（质量值）+33之后的值对应的字符
033!
134"
1043+
20535
3063?
4073I
4174J

为什么偏要+33，而不是+11啊！+22啊之类的？

你可以看看上面的ASCII码表，因为在0~32这个范围之内有很多奇形怪状的字符，同时还包含什么制表符、换行符、回车符、空格、响铃、退格等这类的看不见的字符。而很多时候质量值需要一是肉眼能方便查看；二是要方便程序读取。显然在0-32这个范围内不适合取。而从33开始直到127均为可见的字符。那从33开始取值不就的了，怎么又来个什么+64之类的？

谷歌了一下，暂时没有找到究竟+64的原因，但是在一篇博客《质量值体系 Phred33 和 Phred 64 的由来 及其在质量控制中的实际影响 - Part 1》中作者自己脑补的两者之间的关联。我觉得下面这句话可能可以说明一下作者脑补的内容

一流的企业做标准，二流的企业做品牌，三流的企业做产品

就是当时测序公司有多家，互相之间相互竞争，各家公司是竞争对手，怎么可能你用什么我也用什么。所以我也要想出一定的标准，哈哈！后来的Solexa就不采用+33的方式（这个方式就是Phred+33的编码方式当时Sanger公司所采用的），而使用Solexa+64，其实+64也不是就是随便来的，可以再看一下上面的ASCII码对照表。

Quality encoding: @ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefgh
                  |         |         |         |         |
   Quality score: 0........10........20........30........40

在质量值0对应的@之后，就是从英文大写字母A开始了，也就是质量值1-26对应编码字符A-Z了。这样人肉眼看的话其实质量值是相当直观了。不过在质量值到达33之后，又出现了小写字母......这个个人感觉就不如Phred+33那样的呢。

到后来2006年Solexa公司被Illumina公式收购，于是又有了Illumina 1.3+ Phred+64（这里是我自己脑补的- . -）。于是就有了Sanger公司的Phred+33的方式以及Illumina的Phred+64的方式。于是两家公司互相不让着谁。这可苦了当时编写生信数据分析的人，后来没办法只有在程序中加上自动判断的程序。

可是墨菲他老人家不高兴了（这里我不太了解墨菲和这个事件之间的联系），所以在 2011 年, Illumina 公司表示他们又要改成 Phred33 体系了 (Upcoming changes in CASAVA)。这样来来回回还是回到了Phred+33的体系。

期间有意思的是，当时三巨头的另一家测序仪公司 454 Life Sciences (后被 Roche 收购) 就更绝, 人家从碱基开始就不用 ACTG 表示, 直接整了个 ColorSpace 体系出来, 根本不和你们玩，当然后来大家也不跟 454 玩了, 最后他也就没得玩了 。

这个是在网上别人得到的 ColorSpace数据

@SRR2967009.1 100_1000_1168_F3
T10011023211201220121202030102221012302121010131001
+2@@@@>@?@@@@<@@//;@@/@9?@8@=@@@6;6@66;<@6@67?2?;/@@SRR2967009.2 100_1000_1211_F3
T20132312201120021312220200023110220113100012321011
+
@@@@@@@@@<@@@@@@@@@@@@@@@@@@@@@@?@@@@/?@@@@@@@@<?@

序列是用double-encoded data或者称colorspace的模式表示的。

按照手册的说明

AA, CC, GG, TT : 0
AC, CA, GT, TG : 1
AG, CT, GA, TC : 2
AT, CG, GC, TA : 3

Therefore five base long sequence AACTA will be represented as 0123. The encoding is AA=0, AC=1, CT=2, and TA=3.

额，感觉有些诡异。

2.2.0.5 质量值的区别方法

  1. 肉眼查看
    通过Phred+33与Phred+64的区别，可以知道，[0~9]的为Phred33特有的质量字符，小写字母[a~z]的为Phred64特有的质量字符。

  2. 程序判断

这个shell代码不知道出处是哪里，反正放在这里吧，如果有侵权，那我就删除了。

fqtype () {
        less $1 | head -n 999 | awk '{if(NR%4==0) printf("%s",$0);}' \
        | od -A n -t u1 -v \
        | awk 'BEGIN{min=100;max=0;}\
        {for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END\
        {if(max<=74 && min<59) print "Phred+33"; \
        else if(max>73 && min>=64) print "Phred+64"; \
        else if(min>=59 && min<64 && max>73) print "Solexa+64"; \
        else print "Unknown score encoding"; \
        print "( " min ", " max, ")";}'}

仿照上面的shell代码，我改成了Perl One-liners的版本。一般只需要取一部分序列进行计算就可以了，计算这个值还是挺快的。使用真实的fastq数据同样也很快，只需要接近0.05秒的时间

# 这里实际上只取20行不够，真正的最好取上千条比较好(4000行)# 因为有的时候测序质量很高，采样太少根本无法得到真正的上下限的值# 这里是因为123.fastq.gz文件本身只有80行，所以只取了80行# 为了方便使用，这里可以提前赋值一个变量，以后改行数的时候会比较方便test_line=80
zcat 123.fastq.gz | head -n $test_line | perl -n -e '
  BEGIN{
    # 首先初始化两个用来判断质量值上下限的两个变量
    $max = 0;    # 将最大值设置得精良小
    $min = 1000; # 将最小值设置得尽量大
  }
  # 质量行在4的倍数行
  if($. % 4 == 0){
    # 去除回车符换行符
    s/\r?\n$//;
    my $word;
    # 每次“吞”一个字符进行ASCII码值转换
    while($word = chop($_)){
      my $ascll_num = ord($word);
      # 不断扩大min与max的范围，得到更加真实的上下限值
      if($ascll_num > $max){$max = $ascll_num;}
      if($ascll_num < $min){$min = $ascll_num;}
    }
  }
  END{
    # 如果质量最大值低于75，最小值小于59，那么认为就是Phred+33类型
    if($max < 75 and $min < 59){
      print "Phred+33!\n";
    }elsif($max > 73 and $min > 63){
    # 如果质量最大值高于73，最小值高于63，那么认为就是Phred+64类型
      print "Phred+64!\n";
    }elsif($min > 53 and $min < 64 and $max > 73){
    # 如果最小值在58和64之间，最大值大于73，那就是Solexa+64类型
      print "Solexa+64!\n";
    }else{
    # 否则不知道类型
      print "unkown!\n";
    }
    # 打印出上下限信息
    print "Max:$max\n";
    print "Min:$min\n";
  }
'

输出为

Phred+33!Max:74Min:35

但是这里你觉得这样的步骤可不可以分解呢？

# [说明]：# 这里不是炫技之类的，就是多使用其他linux工具进行搭配# 减少perl中的那些判断的代码，让单行程序更加“单行”一些# 让功能更加专一，便于代码重用# zcat 负责解压# head 负责提取前面几行# awk  负责排除其他行# perl 负责得到上下限的值# perl 负责对这个值进行判断并且得到质量类型test_line=80
zcat 123.fastq.gz | head -n $test_line | awk 'NR%4==2{print}' | perl -n -e '
  BEGIN{
    # 首先初始化两个用来判断质量值上下限的两个变量
    $max = 0;    # 将最大值设置得尽量小
    $min = 1000; # 将最小值设置得尽量大
  }
  # 去除回车符换行符
  s/\r?\n$//;
  my $word;
  # 每次“吞”一个字符进行ASCII码值转换
  while($word = chop($_)){
    my $ascll_num = ord($word);
    # 将上下限的值赋值给刚才的两个变量
    if($ascll_num > $max){$max = $ascll_num;}
    if($ascll_num < $min){$min = $ascll_num;}
  }
  END{
    print "$max:$min\n";
  }
' | perl -n -a -F: -e '
  BEGIN{
    sub print_type {
      my ($max,$min)=@_;
      if($max < 75 and $min < 59){
        return "Phred+33!\n";
      }elsif($max > 73 and $min > 63){
        return "Phred+64!\n";
      }elsif($min > 53 and $min < 64 and $max > 73){
        return "Solexa+64!\n";
      }else{
        return "unkown!\n";
      }
    }
  }
  my ($max,$min) = @F;
  print "Max:$max\n";
  print "Min:$min\n";
'

2.3 fastqc相关统计量[选看]

2.3.1 序列长度的分布

这个统计序列长度就比统计fasta序列的长度简单一些啦！

# 首先使用awk排除其他行zcat 123.fastq.gz | awk 'NR%4==2{print}' | perl -n -e '
  s/\r?\n$//;
  my $seq_len = length($_);
  # 把序列长度以及为此长度的信息存入哈希中
  $hash{$seq_len}++;
  END{
    # 打印出不同长度以及为此长度的序列的条数
    print "length\tcount\n";
    for my $len (sort {$a <=> $b} keys %hash){
      print "$len\t$hash{$len}\n";
    }
  }
'

输出为

length  count
89      1
90      18
91      1

也可以重定向到文件中，拿这个文件可以使用R来画图。

同样的你也可以写成

zcat 123.fastq.gz | perl -n -e '
  if($. % 4 == 2){
    s/\r?\n$//;
    my $seq_len = length($_);
    # 把序列长度以及为此长度的信息存入哈希中
    $hash{$seq_len}++;
  }
  END{
    # 打印出不同长度以及为此长度的序列的条数
    print "length\tcount\n";
    # 按照长度进行排序
    for my $len (sort {$a <=> $b} keys %hash){
      print "$len\t$hash{$len}\n";
    }
  }
'

实际上上面的方法比下面的要慢一些，但是上面的程序将perl所承担的任务量减少了，所以代码量看起来少了一些。有的时候不能兼得，每项需要均衡一下。

放松时间

|_・) |ω・） |･ω･｀) (╯ε╰) ^.^ (๑￫ܫ￩) _ (╯ε╰) ^.^ (๑￫ܫ￩) _ (╯ε╰) ^.^ (๑￫ܫ￩) _ ﾍ(･_|

最近听了一下重阳世界观-狂点技能树(1)，和大家分享一下，里面说到：

在兵马俑坑，出土最多的青铜兵器是箭头，考古人员发现了4万多支箭头，是三棱型的。制作的非常规范，箭头底部宽度的误差为0.83毫米，而且金属配比基本上是相同的。

为什么选择这种三棱形状的箭头呢？

三棱箭头拥有三个锋利的棱角，在击中目标的瞬间，棱的锋刃处就会形成切割力，箭头就能够穿透铠甲，直达人体。同时这三个面是呈流线型的，据说这个流线型的轮廓和子弹的外形都一致。按照现在话来说，有一个好的空气动力学的特性。除了三棱型的箭矢，还有什么狼牙箭——带倒钩和翼面的。但翼面容易受风的影响，使箭头偏离目标。

武器有很多个维度，如果只追求单一的维度，那么这个武器不会太好，就说这样一个小小的箭头啊！制作工艺简单，标准化，那样才能大量生产；如果太复杂，成本太高，不能大量装备；另外如果只盲目追求杀伤力，打不准那也不行。

三棱箭头.jpg


各种各样的箭头.jpg

2.3.2 每个位置上质量值的分布

其实在fastQC软件之中就对这个做了一个很好的诠释，这里我用Perl代码+R代码尝试实现这种每个位置上的碱基质量分布

zcat 123.fastq.gz | perl -n -e '
  BEGIN{
    # 先声明一个质量值转换的子程序
    sub convert_quality_to_score{
      my $ascll = shift;
      # 这里我默认使用Phred+33
      return ord($ascll) - 33;
    }
  }
  # 只需要每个单位的第四行
  if($. % 4 == 0){
    s/\r?\n//;
    my $n = 0;
    while($_){
      $n++;
      my $quality_score = substr($_,0,1,"");
      # 将碱基推到对应的位置的哈西中
      push @{$hash{$n}},convert_quality_to_score($quality_score);
    }
  }
  END{
    my @list;
    # 设置列表内插之后打印元素的分隔符
    $" = ",";
    for my $base_site (sort {$a <=> $b} keys %hash){
      # 按照R语言的方式写入信息
      print "c_${base_site} <- c(@{$hash{$base_site}})\n";
      push @list,"c_${base_site}";
    }
    # 写入执行箱线图的boxplot
    print "boxplot(@list)\n";
  }
' > statistic_base_quality.R# 打开命令行的R# 这一步windows需要安装一下R，并且设置环境变量R# 执行R脚本source("statistic_base_quality.R")



作者：白菜代码小推车
链接：https://www.jianshu.com/p/6066d471bcbe


相关标签
Java
0人推荐
随时随地看视频
慕课网APP