通过简单数据熟悉Linux下生物信息学各种操作4

原地址
几点说明
1.非简单翻译，所有代码均可运行，为了辅助理解，基本每步代码都有结果，需要比较的进行了整合
2.原文中的软件都下载最新版本
3.原文中有少量代码是错误的，这里进行了修正
4.对于需要的一些知识背景，在这里进行了注释或链接到他人博客

一共4部分
通过简单数据熟悉Linux下生物信息学各种操作1
通过简单数据熟悉Linux下生物信息学各种操作2
通过简单数据熟悉Linux下生物信息学各种操作3
通过简单数据熟悉Linux下生物信息学各种操作4

20 Pileup和Coverage

上面已经得到的mutation files
需要用到前面的内容

20.1 安装bcftools suite以便处理vcf文件（variant call formats）

cd src
curl -OL http://sourceforge.net/projects/samtools/files/samtools/1.9/bcftools-1.9.tar.bz2
tar jxvf bcftools-1.9.tar.bz2
cd bcftools-1.9
make
ln -s ~/src/bcftools-1.9/bcftools ~/bin

每个run的平均coverage是什么
默认，samtools depth只输出非0 coverage的区域

samtools depth bwa.bam | awk '{ sum += $3 } END { print sum/NR } '

73.8133

以下内容部分和前面重复

20.2 计算genome大小

samtools view -h bow.bam |head| samtools view -h bow.bam |head|grep @SQ

@SQ SN:NC_002549    LN:18959

samtools depth bwa.bam | awk '{ sum += $3 } END { print sum/18959 } '

结果仍然是

73.7938

20.3 分别看有何无参考基因组的pileups

前已述及，不再展示图

samtools mpileup -Q 0 bwa.bam | more
samtools mpileup -f ~/refs/852/NC.fa -Q 0 bwa.bam | more

20.4samtools tview（samtools文本基因组浏览器）

samtools tview bwa.bam

1         11        21        31        41        51        61        71        
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    CACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
    CACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTG
      CACAAAAAGAAAGAAGAATTTTTAGGATCTTTAGTGTGCGAATAACTATGAGGAATATTAATAATTTACCTCTC
          AAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
           AAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
            AAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
              GAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAAGTATGAGGAAGATTAATAATTTTCCTCTC
                 AGAAGAATTTTTAGGATCTTTTGTGTTCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
                   AAGAATTTTTAGGATCTTTTGTGTGGGAATAACTATGAGGAAGATTCAAAATTTTCCTCTC
                    AGAATTTTTAGTATCTTTTGAGTGCGACTAACTATGAGGAAGATTAATAATTTTCCTCTC
                      AATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCAC
                      AATTTTTAGGATCTTTTGTGTGCGAATCACTATGAGGAAGATTAATAATTTTCCTCTC
                       ATTATTAGGATCTTTTGTGTGCGAATAACTATGAGGACGATTAATAATTTTCCTCTC
                        TTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATACTTTTCCTCTC
                          TTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
                           TTAGGATCTTTTGTGTGCGAATAACTATGATGAAGATTAATAATTTTCCTCTC
                             AGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
                                   TATTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
                                   TTTTGTGTGCGAATAAGTATGAGGAAGATTAATAATTTTCCTCTC
                                    TTTGTGTGCGAATAACTATGAGGAAGATTAATCATTTTCCTCTC
                                     ATGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
                                         GTGCGAATAAGTATGCGGCAGATTAATAATTTTCCTCTC
                                                TAACTATGAGGAAGATTAATAATTTTCCTCTC
                                                TAACTATGAGGAAGATTAATAATTTTCCTCTC
                                                   CTATGAGGAAGATTAATAATTTTGCTCTC
                                                     ATGAGGAAGATTAATAATTTTCCTCTC
                                                      TGCGGAAGATTAATAATTTTCCTCTC
                                                        AGGAAGATAAATAATTTTCCTCTC
                                                            AGATTAATAATTTTCCTCTC
                                                            AGATTAATAATTTTCCTCTC
                                                               TTAATAATTTTCCTCTC
                                                                     ATTTTCCTCTC
                                                                      TTTTCCTCTC
                                                                      TTTTCCTCTC
                                                                       TTTCCTCTC
                                                                            TCTC
                                                                             CTG
                                                                               C

20.5 对整个genome生成vcf

samtools mpileup -Q 0 -f ~/refs/852/NC.fa -uv bwa.bam | more

部分结果如下

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  liver
NC_002549       5       .       C       <*>     0       .       DP=1;I16=1,0,0,0,17,289,0,0,60,3600,0,0,0,0,0,0;QS=1,0;MQ0F=0   PL      0,3,17
NC_002549       6       .       A       <*>     0       .       DP=1;I16=1,0,0,0,17,289,0,0,60,3600,0,0,1,1,0,0;QS=1,0;MQ0F=0   PL      0,3,17
NC_002549       7       .       C       <*>     0       .       DP=2;I16=2,0,0,0,34,578,0,0,120,7200,0,0,2,4,0,0;QS=1,0;MQ0F=0  PL      0,6,31
NC_002549       8       .       A       <*>     0       .       DP=2;I16=2,0,0,0,34,578,0,0,120,7200,0,0,4,10,0,0;QS=1,0;MQ0F=0 PL      0,6,31
NC_002549       9       .       C       <*>     0       .       DP=2;I16=2,0,0,0,34,578,0,0,120,7200,0,0,6,20,0,0;QS=1,0;MQ0F=0 PL      0,6,31
NC_002549       10      .       A       <*>     0       .       DP=2;I16=2,0,0,0,34,578,0,0,120,7200,0,0,8,34,0,0;QS=1,0;MQ0F=0 PL      0,6,31
NC_002549       11      .       A       <*>     0       .       DP=3;I16=3,0,0,0,51,867,0,0,180,10800,0,0,10,52,0,0;QS=1,0;MQ0F=0       PL      0,9,42
NC_002549       12      .       A       <*>     0       .       DP=4;I16=4,0,0,0,68,1156,0,0,240,14400,0,0,13,75,0,0;QS=1,0;MQ0F=0      PL      0,12,50
NC_002549       13      .       A       <*>     0       .       DP=5;I16=5,0,0,0,85,1445,0,0,300,18000,0,0,17,105,0,0;QS=1,0;MQ0F=0     PL      0,15,57
NC_002549       14      .       A       <*>     0       .       DP=5;I16=5,0,0,0,85,1445,0,0,300,18000,0,0,22,144,0,0;QS=1,0;MQ0F=0     PL      0,15,57
NC_002549       15      .       G       <*>     0       .       DP=6;I16=6,0,0,0,102,1734,0,0,360,21600,0,0,27,193,0,0;QS=1,0;MQ0F=0    PL      0,18,62
NC_002549       16      .       A       <*>     0       .       DP=6;I16=6,0,0,0,102,1734,0,0,360,21600,0,0,33,253,0,0;QS=1,0;MQ0F=0    PL      0,18,62
NC_002549       17      .       A       <*>     0       .       DP=6;I16=6,0,0,0,102,1734,0,0,360,21600,0,0,39,325,0,0;QS=1,0;MQ0F=0    PL      0,18,62
NC_002549       18      .       A       <*>     0       .       DP=7;I16=7,0,0,0,119,2023,0,0,420,25200,0,0,45,409,0,0;QS=1,0;MQ0F=0    PL      0,21,66
NC_002549       19      .       G       <*>     0       .       DP=7;I16=7,0,0,0,119,2023,0,0,420,25200,0,0,52,506,0,0;QS=1,0;MQ0F=0    PL      0,21,66
NC_002549       20      .       A       <*>     0       .       DP=8;I16=8,0,0,0,136,2312,0,0,480,28800,0,0,59,617,0,0;QS=1,0;MQ0F=0    PL      0,24,69
NC_002549       21      .       A       <*>     0       .       DP=9;I16=9,0,0,0,153,2601,0,0,540,32400,0,0,67,743,0,0;QS=1,0;MQ0F=0    PL      0,27,72

20.6 生成genotypes然后call variants

samtools mpileup -Q 0 -ugf ~/refs/852/NC.fa bwa.bam | bcftools call -vc > samtools.vcf

查看snp calls，知道snp calling的原理
这里原作者写了一个python脚本名字为snpcaller.py.但我还没找到
代码先放这里，但不影响后续分析。做个

记号

cat bwa.pileup | python snpcaller.py > diy.vcf

cat samtools.vcf |grep -v "##"|cut -f 1-5|head

#CHROM  POS ID  REF ALT
NC_002549   425 .   G   T,A
NC_002549   507 .   G   T
NC_002549   2846    .   a   aT
NC_002549   2847    .   g   gT
NC_002549   7894    .   A   C
NC_002549   9569    .   G   A
NC_002549   12101   .   T   A
NC_002549   12957   .   G   A
NC_002549   14086   .   A   G

20.7安装freebays以call variants

首先安装Freebays，mac需要另外一个cmake才可以对freebays进行编译,所以先安装cmake

brew install cmakebrew install cmake

安装freebays
关于freebays，可以看下下面的文章
In-depth-NGS-Data-Analysis-Course
Calling variants with freebayes

cd ~/src
git clone --recursive git://github.com/ekg/freebayes.git
cd freebayes
make

注意原文中说，这个地方有个bug，make第一次会报错，然后再make一次就OK。我这里还是显示如下错误

fatal error: 'lzma.h' file not found

通过如下方式解决

brew install xz

做链接

ln -sf ~/src/freebayes/bin/freebayes ~/bin

20.8用安装的FreeBayes call snps

freebayes -f ~/refs/852/NC.fa bwa.bam > freebayes.vcf

对结果进行可视化。关于vcf格式请参考这篇文章VCF格式

Lecture 22变异效应预测

安装bioawk

cd ~/src
git clone https://github.com/lh3/bioawk.git
cd bioawk
make
ln -sf ~/src/bioawk/bioawk ~/bincd ~/src
git clone https://github.com/lh3/bioawk.git
cd bioawk
make
ln -sf ~/src/bioawk/bioawk ~/bin

bioawk的manual在这里

man ~/src/bioawk/awk.1

bioawk能辨识生物信息学类型。

cat r1.fq | bioawk -c fastx ' { print $seq } ' | head -1
cat mutations.gff | bioawk -c gff ' { print $seqname, $start } ' | head -1
cat mutations.gff | bioawk -c gff ' { print $seqname, $end-$start+1 } ' | head -1

下载snpEff,注意原文中下载地址有错

curl -OL https://sourceforge.net/projects/snpeff/files/snpEff_latest_core.zip/download

可以向trimmomatic和readseq一样设置snpEff，现在用alias的方式设置

alias snpEff='java -jar ~/src/snpEff/snpEff.jar'
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 222,681评论 6赞 517
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 95,205评论 3赞 399
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 169,421评论 0赞 362
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 60,114评论 1赞 300
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 69,116评论 6赞 398
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 52,713评论 1赞 312
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 41,170评论 3赞 422
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 40,116评论 0赞 277
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 46,651评论 1赞 320
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 38,714评论 3赞 342
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,865评论 1赞 353
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 36,527评论 5赞 351
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 42,211评论 3赞 336
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,699评论 0赞 25
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,814评论 1赞 274
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 49,299评论 3赞 379
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 45,870评论 2赞 361