前往小程序，Get更优阅读体验！
立即前往
文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
首页
学习
活动
专区
工具
TVP最新优惠活动
返回腾讯云官网
社区首页 >专栏 >通过简单数据熟悉Linux下生物信息学各种操作4

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

作者头像
Y大宽
发布于 2019-07-03 17:15:16
8020
发布于 2019-07-03 17:15:16
举报
文章被收录于专栏：Y大宽

原地址

几点说明

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）

代码语言：javascript
复制
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的区域

代码语言：javascript
复制
samtools depth bwa.bam | awk '{ sum += $3 } END { print sum/NR } '
代码语言：javascript
复制
73.8133

以下内容部分和前面重复

20.2 计算genome大小

代码语言：javascript
复制
samtools view -h bow.bam |head| samtools view -h bow.bam |head|grep @SQ
代码语言：javascript
复制
@SQ SN:NC_002549    LN:18959
代码语言：javascript
复制
samtools depth bwa.bam | awk '{ sum += $3 } END { print sum/18959 } '

结果仍然是

代码语言：javascript
复制
73.7938

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

前已述及，不再展示图

代码语言：javascript
复制
samtools mpileup -Q 0 bwa.bam | more
samtools mpileup -f ~/refs/852/NC.fa -Q 0 bwa.bam | more

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

代码语言：javascript
复制
samtools tview bwa.bam
代码语言：javascript
复制
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

代码语言：javascript
复制
samtools mpileup -Q 0 -f ~/refs/852/NC.fa -uv bwa.bam | more

部分结果如下

代码语言：javascript
复制
#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

代码语言：javascript
复制
samtools mpileup -Q 0 -ugf ~/refs/852/NC.fa bwa.bam | bcftools call -vc > samtools.vcf

查看snp calls，知道snp calling的原理

这里原作者写了一个python脚本名字为snpcaller.py.但我还没找到

代码先放这里，但不影响后续分析。做个

记号

代码语言：javascript
复制
cat bwa.pileup | python snpcaller.py > diy.vcf
代码语言：javascript
复制
cat samtools.vcf |grep -v "##"|cut -f 1-5|head
代码语言：javascript
复制
#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
本文参与 腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。
原始发表：2019.07.03 ，如有侵权请联系 cloudcommunity@tencent.com 删除
linux

本文分享自 作者个人站点/博客 前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ，欢迎热爱写作的你一起参与！

linux
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
LV.
文章
0
获赞
0
目录
  • 20 Pileup和Coverage
    • 20.1 安装bcftools suite以便处理vcf文件（variant call formats）
      • 20.2 计算genome大小
        • 20.3 分别看有何无参考基因组的pileups
          • 20.4samtools tview（samtools文本基因组浏览器）
            • 20.5 对整个genome生成vcf
              • 20.6 生成genotypes然后call variants
              • 记号
              领券

              腾讯云开发者

              扫码关注腾讯云开发者

              扫码关注腾讯云开发者

              领取腾讯云代金券

              热门产品

              热门推荐

              更多推荐

              Copyright © 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

              深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

              腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

              问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

              Copyright © 2013 - 2024 Tencent Cloud.

              All Rights Reserved. 腾讯云 版权所有

              登录 后参与评论