七 17

转-基因突变种类大全

Posted on by

突变（Mutation, 即基因突变）:在生物学上的含义，是指细胞中的遗传基因（通常指存在于细胞核中的脱氧核糖核酸）发生的改变。它包括单个碱基改变所引起的点突变，或多个碱基的缺失、重复和插入。原因可以是细胞分裂时遗传基因的复制发生错误、或受化学物质、辐射或病毒的影响。

以功能分类:

失去功能的突变Loss-of-function mutations

：失去功能的突变是指发生的突变会造成基因完全地失去活性，原因可分成两类。一类是由于基因被删除或是调控基因表现的过程受到影响让基因不表现，另一种则是由于基因本身受到影响，使得基因的产物蛋白质失去功能。又称剔除突变（null mutations）或是敲除突变（knockout mutations）。

次形态突变Time form mutation：此种突变会使基因的表现或是基因产物的活性减弱，但不会消失。

超形态突变hypermorphic mutations：此种突变与次形态突变相反，会使基因的表现加强。

获得功能的突变gain-of-function mutation：获得功能的突变是指发生的突变让原本应该是不表现的基因产生活性，进而影响细胞功能，这样的突变多半需要染色体程度的突变较有可能产生，而最常发生获得功能的突变就是癌细胞。

以突变机理分类:

  1. 点突变point mutation:DNA序列中涉及单个核苷酸或碱基的变化称为点突变。 通常有两种情况：一是一种碱基或核苷酸被另一种碱基或核苷酸所替换；二是一个碱基的插入或缺失。

                   （1）沉默突变silent mutation

当点突变发生在基因及其调控序列之外，或使基因序列内一种密码子变成编码同一种氨基酸的另一种同义密码子时，不会改变生物个体的基因产物，因而不引起性状变异。不引起生物性状变异的突变称为沉默突变。

                   （2）错义突变missense mutation

指由于某个碱基对的改变，使编码一种氨基酸的密码子变成编码另外一种氨基酸的密码子，结果使构成蛋白质的数百上千个氨基酸中有一个氨基酸发生变化。（实例：镰刀形细胞贫血症）

                   （3）移码突变frameshift mutation

指在DNA链上，有时一个或几个非3的整数倍的碱基的插入或缺失，往往产生比碱基替换突变更严重的后果。 这种插入或缺失突变会造成阅读框的改变，翻译过程中其下游的三联密码子都被错读，产生完全错误的肽链或肽链合成提前终止。这种插入或缺失突变又称为移码突变。

                   （4）无义突变nonsense mutation

是指当点突变使一个编码氨基酸的密码子变成终止子时，则蛋白质合成进行到该突变位点时会提前终止，结果产生一个较短的多肽链或较小的蛋白质。

  1. 大突变

大突变是可能涉及整个基因以至多个基因的一长段DNA序列的改变，大突变常常导致染色体畸变。

（1）缺失:指DNA分子丢失一段碱基序列。（染色体缺失）(Deletion)

（2）插入:指DNA分子的正常序列中插入一段DNA序列。(Insertion②)

（3）重排:重排包括某段DNA序列的重复（duplication），倒位（inversion），易位（translocation）等。

 

七 17

转-R语言内存管理

Posted on by

R中的对象（比如矩阵）在内存中存于两种不同的地方:

第一种是堆内存（heap），其基本单元是“Vcells”，每个大小为8字节，新来一个对象就会申请一块空间，把值全部存在这里，和C里面的堆内存很像;

第二种是地址对（cons cells），主要用来存储地址信息，最小单元一般在32位系统中是28字节、64位系统中是56字节。

1、ls()来查看当前所有对象名，对于每一个对象，可以通过object.size(x)来查看其占用内存的大小。

如果是因为当前对象占用内存过多，那么可以通过处理对象来获取更大的可用内存。一个很有用的方法是改变对象的存储模式，通过storage.mode（x）可以看到某个对象的存储模式，比如某个矩阵默认就是“double”的，如果这个矩阵的数值都是整数甚至0-1，完全没必要使用double来占用空间，可以使用storage.mode(x） <- "integer"将其改为整数型，可以看到该对象的大小会变为原来的一半。

2、object.size()看每个变量占多大内存。

3、memory.size()查看现在的work space的内存使用

4、memory.limit()查看系统规定的内存使用上限。如果现在的内存上限不够用，可以通过memory.limit(newLimit)更改到一个新的上限。注意，在32位的R中，封顶上限为4G，无法在一个程序上使用超过4G （数位上限）。这种时候，可以考虑使用64位的版本。

 

对于一些很大的但无用的中间变量，养成清理的习惯：

可以使用rm(object)删除变量，但是记住，rm后记得使用gc()做Garbage collection，否则内存是不会自动释放的，相当于你没做rm.

七 16

模拟测序lambda_virus基因组

Posted on by
lambda_virus基因组文件是bowtie软件自带的测试数据，共48502个bp，首先我用脚本模拟出它的全打断文件！

perl -alne '{next if /^>/;$a.=$_;}END{$len=length $a;print substr($a.$a,$_,120) foreach 0..$len}' lambda_virus.fa >lamb_virus.120bp

长度均为120bp的片段。

我测序的策略是CTAG碱基重复30次，共加入120个碱基。

对每个120bp片段来说，如果遇到互补碱基就加上，直到120个碱基加完，这样如果比较巧合的话，会有部分碱基能全部加满120bp的，但是如果每个120bp片段的ATCG分布均匀，那么就都应该30bp碱基能被加上。

image001

[perl]
while (<>) {

$seq=$_;$sum=0;

foreach $i (0..120){

$str=substr($seq,$i,2);

if ($str eq "GG"| $str eq "CC"| $str eq "AA"| $str eq "TT"){$sum+=4;}

elsif ($str eq "GT"| $str eq "CG"| $str eq "AC"| $str eq "TA"){$sum+=3;}

elsif ($str eq "GA"|$str eq "CT"| $str eq "AG"| $str eq "TC"){$sum+=2;}

else{$sum+=1;};

#print "$sum\n";

if ($sum>120){print "$i\n";last;}

}

}

[/perl]

perl length.pl lambda_virus.120bp >length.txt

得到结果如下：

 

Length 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
Count 2 19 34 110 204 432 878 1495 2237 3202 4343 5179 5697 5429 4865 4214
Length 52 53 54 55 56 57 58 59 60 61 62 63 64
Count 3249 2499 1735 1090 657 396 228 141 90 48 18 9 3

右表可以看出，大部分测序得到碱基长度都集中在46bp到51bp之间

画出箱线图如下

image003

画出条形图如下：

image005

 

然后我模拟了一个6000bp的基因组，做同样的模拟测序看看评价测序长度分布情况：

Length 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Count 9 22 96 207 322 382 479 671 770 714 706 546 424 232 182 100 52 30 14 21
Length 59 60 61                                  
Count 15 5 2                                  

可以看出这次的测序片段集中在45到51bp

七 15

perl操作excel表格

Posted on by

perl这个语言已经过时很久了，所以它的模块支持性能不是很好，暂时我只看到了对excel2003格式的表格的读取写入操作。

以下是我参考Spreadsheet::ParseExcel这个模块写的一个把excel表格转为csv的小程序，大家也可以自己搜索该模块的说明文档，这样学的更快一点！

[perl]

#!/usr/bin/perl -w
# For each tab (worksheet) in a file (workbook),
# spit out columns separated by ",",
# and rows separated by c/r.

use Spreadsheet::ParseExcel;
use strict;
use utf8;
use Encode::Locale qw($ENCODING_LOCALE_FS);
use Encode;

my $filename ="test.xls";#输入需要解析的excel文件名，必须是03版本的
my $e = new Spreadsheet::ParseExcel; #新建一个excel表格操作器
my $eBook = $e->Parse($filename);    #用表格操作器来解析我们的文件
my $sheets = $eBook->{SheetCount};   #得到该文件中sheet总数
my ($eSheet, $sheetName);

foreach my $sheet (0 .. $sheets - 1) {
$eSheet = $eBook->{Worksheet}[$sheet];
$sheetName = $eSheet->{Name};
my $f1 = encode(locale_fs => $sheetName); #每次操作中文我都很纠结，还得各种转码
open FH_out ,">$f1.csv" or die "error open ";
next unless (exists ($eSheet->{MaxRow}) and (exists ($eSheet->{MaxCol})));
foreach my $row ($eSheet->{MinRow} .. $eSheet->{MaxRow}) {
foreach my $column ($eSheet->{MinCol} .. $eSheet->{MaxCol}) {
if (defined $eSheet->{Cells}[$row][$column])
{
print FH_out $eSheet->{Cells}[$row][$column]->Value . ",";
} else {
print FH_out ",";
}
}
print FH_out "\n";
}
close FH_out;

}
exit;

[/perl]

 

七 13

使用Seq2HLA进行HLA分型

Posted on by

基于高通量测序数据进行HLA分型的软件挺多的，比较老的有三个，作者分别是Boegel et al.Kim et al.Major et al.，然后他们都被OptiType这个软件的作者被批评了，我这里先介绍Kim et al的seq2HLA使用方法，以下是它的一些链接。

功能概述

seq2HLA is a computational tool to determine Human Leukocyte Antigen (HLA) directly from existing and future short RNA-Seq reads. It takes standard RNA-Seq sequence reads in fastq format as input, uses a bowtie index comprising known HLA alleles and outputs the most likely HLA class I and class II types, a p-value for each call, and the expression of each class.

软件简介

Type of tool     Program

Nature of tool          Standalone

Operating system   Unix/Linux, Mac OS X

Language        Python, R

Article     (Boegel et al., 2013) HLA typing from RNA-Seq sequence reads. Genome medicine.

PubMed http://www.ncbi.nlm.nih.gov/pubmed/23259685

URL          https://bitbucket.org/sebastian_boegel/seq2hla

源代码，下载并安装

https://bitbucket.org/sebastian_boegel/seq2hla/src

http://tron-mainz.de/tron-facilities/computational-medicine/seq2hla/

第一版是这样的

image001

第二版是这样的

image002

只有第二版才支持gz压缩包格式的fastq，而且不需要指定length了

其中reference文件夹下面的是发布这个软件的团体已经制备好来的HLA库文件

image003

下载即可使用，前提是你的系统其它环境都OK

用法：

python seq2HLA.py -1 <readfile1> -2 <readfile2> -r "<runname>" [-p <int>]* [-3 <int>]**

image004

很简单，-1和-2指定我们的双端测序数据即可，可以是压缩包格式的（自动调用gzip），-r的输出目录，会输出7个文件，需要一个个解读，-p指定线程数给bowtie用的，-3是指定需要trim几个低质量碱基。

但是运行这个软件的要求非常多，需要安装好python和R，而且还有版本限制，需要安装好biopython而且还必须是双端测序，而且当前文件夹下面的reference文件夹下面必须有参考基因组的bowtie索引，而且系统必须安装好了bowtie，还需要在快捷方式里面！

我这里用的是第二版的

image006

所以，我用的也是第二版改进的命令。非常好用,我这里用的是一个外显子测序数据，是hiseq2500测的PE100

python seq2HLA.py -1 ../../6-exon/PC3-1.read1_Clean.fastq.gz -2 ../../6-exon/PC3-1.read2_Clean.fastq.gz -r PC3

貌似输出文件太多了一点

#Output:#The results are output to stdout and to textfiles. Most important are:

#i) <prefix>-ClassI.HLAgenotype2digits => 2 digit result of Class I

#ii) <prefix>-ClassII.HLAgenotype2digits => 2 digit result of Class II

#iii) <prefix>-ClassI.HLAgenotype4digits => 4 digit result of Class I

#iv) <prefix>-ClassII.HLAgenotype4digits => 4 digit result of Class II

#v) <prefix>.ambiguity => reports typing ambuigities (more than one solution for an allele possible)

#vi) <prefix>-ClassI.expression => expression of Class I alleles

#vii) <prefix>-ClassII.expression => expression of Class II alleles

根据文献，我简单看了一下，文件的确好复杂，不过我们只需要看输出日志即可

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

A       A*68   7.287148e-05   A*24   0.03680272

B       B*52   0.1717737       B*53   0.3952319

C       C*12   0.03009331     hoz("C*14")     0.6783964

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassI.bowtielog

A: 7.93 RPKM

C: 9.75 RPKM

B: 8.35 RPKM

The digital haplotype is written into BC1-1/BC1-1-ClassI.digitalhaplotype3

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!A     A*68:01 7.287148e-05   A*24:02 0.03680272

!B     B*52:01 0.1717737       B*53:01'       0.6542288

!C     C*12:02 0.03371717     C*12:02 0.6783964

上面的HLA的class I的数据结果

接下来是class II的数据结果，是不是很简单呀！

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

DQA     DQA1*01 0.1511134       DQA1*02 0

DQB     DQB1*02 0.02321615     DQB1*05 0.42202

DRB     DRB1*15 2.595144e-05   DRB1*07 0.321219

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassII.bowtielog

DQB1: 4.47 RPKM

DRB1: 5.59 RPKM

DQA1: 0.44 RPKM

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!DQA   DQA1*01:02'     0.1511134       DQA1*02:01     0.0

!DQB   DQB1*02:01'     0.02321615     DQB1*05:01     0.42202

!DRB   DRB1*15:02'     2.595144e-05   DRB1*07:01     0.321219

七 06

GATK使用注意事项

Posted on by

GATK这个软件在做snp-calling的时候使用率非常高，因为之前一直是简单粗略的看看snp情况而已，所以没有具体研究它。

这些天做一些外显子项目以找snp为重点，所以想了想还是用起它，报错非常多，调试了好久才成功。

所以记录一些注意事项!

GATK软件本身是受版权保护的，所以需要申请才能下载使用，大家自己去broad institute申请即可。

下载软件就可以直接使用，java软件不需要安装，但是需要你的机器上面有java，当然软件只是个开始，重点是你还得下载很多配套数据，https://software.broadinstitute.org/gatk/download/bundle（ps:这个链接可能会失效，下面的文件，请自己谷歌找到地址哈。），而且这个时候要明确你的参考基因组版本了！！！ b36/b37/hg18/hg19/hg38，记住b37和hg19并不是完全一样的，有些微区别哦！！！
Continue reading →

七 03

毕业生入深户完全指南

Posted on by

第一步：网上个人测评

申请人登录深圳市人力资源保障局官方网站(www.szhrss.gov.cn)，进入“网上办事”--“网上申办”--“深圳市人才引进(毕业生、在职人才引进)测评与申报系统”，注册个人账户，注册成功后通过个人用户登录系统选择 “毕业生接收”，根据系统提示填写个人信息，填报完成后，点击“保存”--点击“按当前填报信息测评”，系统将判断所填报人员是否符合毕业生接收政策并列出符合的政策条款。

也可以直接去测评网址，注册之后填一些信息https://sz12333.gov.cn/rcyj/

Ps:信息填写要真实，填写完了之后等待审核，一般三到五个工作日即可审核完毕，没什么特殊情况都会通过的，如果查看到自己审核通过了就可以进行第二步啦！

 

第二步：上门签订人事代理协议

符合毕业生接收政策的，即可与市人力资源局认可的人力资源代理机构签订个人申办委托办理协议，委托其办理毕业生接收手续。

上门需要带一些必备的资料，如下所示：

 

 

序号 材料名称
2 接收高等院校应届毕业生呈报表（收原件）
3 毕业生推荐表、成绩单（收原件）
4 学历及学位证书（申报时已毕业的验原件,收复印件；申报时未毕业的报到时验原件，不收复印件）
5 身份证（收复印件）
1
户口簿（户籍证明）

 

 

 

以上所有能带原件的都带上，然后所有原件都有复印一份！

代理机构有很多，大家选择自己最方便的， 我去的是深圳市人才交流服务中心（高新区分部） 

这个步骤需要上门，而且还需要排队，很可能需要排队两个到三个小时。还需要交钱，可能是260左右，可以刷卡。

PS：这个步骤因为要请假，所以大家一定要带全资料！！！办理很简单，主要是排队时间太长，办理完了会给你一个回执，你按照回执的提示15个工作日左右即可查看自己是否办理成功！如果成功了就再来一次，拿接收函！！！

 

第三步：用深圳市的接收函在学校拿报到证和户口迁移证

如果你是刚毕业，报到证还没有，那么这一步很简单，委托学校的同学帮忙即可。

如果你已经被开过报到证了（一般是遣返回老家啦），你就需要改派报到证啦！这个改派其实很简单，你需要自己看看你们学校改派流程，委托同学把新的报到证寄给你即可，如果你的档案还在学校就要求学校档案馆把你的档案通过机要传给深圳（15天左右），如果你的档案被遣返回家或者异地，那么你就要打电话去你的档案所在地要求他们帮你把你的档案通过机要传给深圳（15天左右）！

如果你的户口在学校，那么很简单，去你学校弄一个户口迁移证即可。

如果你的户口在老家，那么就麻烦了，还需要农转非什么的，看看你家里人的关系吧！

Ps：用深圳的接收函回学校成功拿到报到证和户口迁移证之后要随时上网查看自己的档案是否到达深圳。

 

第四步：拿介绍信和深圳市入户人员信息卡

这一个步骤不需要排队，在罗湖人才市场，需要身份证，毕业证，学位证，学历验证报告，报到证和户口迁移证原件及复印件各一份，缺一不可！！！

Ps：学历验证报告在学信网即可弄，请保证有效期至少一年以上！！！

第五步：去派出所办理户口身份证

这个需要预约！

这个需要预约！

这个需要预约！

重要的事情说三遍！如果你预约好了，那么你从罗湖人才市场拿到了介绍信和深圳市入户人员信息卡后就可以直接去派出所啦！！！但是如果你没有预约，你就得再等一个星期等到拿到预约时间后才能去派出所办理！

除了需要你在罗湖人才市场拿到了介绍信和深圳市入户人员信息卡，还需要数码照相回执和身份证，以及它们的复印件！！！

Ps：如果你是落户到高新园区派出所，那么你还有个近路，直接去迈瑞警务室也能完成落户流程！

到这里，落户就完成啦！十个工作日之后去派出所拿新的身份证即可！是不是非常简单呀小朋友们！

当然别忘了最后的彩蛋！网上深圳市新引进人才租房补贴系统

https://sz12333.gov.cn/szhr_pubtalent/talent_login.jsp

点击进入，有惊喜。

科未满30周岁、硕士未满35周岁、博士未满40周岁。租房补贴标准为：本科6000/人，硕士9000元/人；博士12000元/人。

 

总结一下：你需要请假三次或者四次，分别是去签人才引进代理协议，再去签人才引进代理协议的地方拿接收函，去罗湖人才市场拿介绍信和入户信息卡，去派出所办理落户及新身份证！

这个流程如果你仔细看了，而且保证按照流程走，当然，以你在各个单位拿到的最新资料为准，记住，各种材料宁可多带，也不能缺，一旦你少带了什么，没有人会跟你讲人情的，一切推倒重来！应该还算是蛮简单的，如果有任何疑问，欢迎咨询我QQ1227278128

六 24

Genomemapper软件使用说明书

Posted on by

 我以前一直以为有了bwa跟bowtie，没什么必要用其它的alignment软件，直到我碰到了高插入删除的helicos三代测序数据，我才发现，这个古董软件genomemapper居然大有用武之地了。

一．下载并且安装该软件

这是最新版本了

Release 0.4.4 2012-10-30 source code including documentation

Wget http://1001genomes.org/data/software/genomemapper/genomemapper_0.4.4/genomemapper-0.4.4.tar.gz

这个软件安装很简单，解压进入目录，make一下即可

image001

看到make完了之后就会多了两个软件，其中一个是用来构建参考基因组索引，一个用来比对的！

二．准备数据

既然是比对软件，那么肯定是一个参考基因组，一个测序的fastq原始文件咯

当然这个软件比较奇葩，它还支持Multi-FASTA, FASTQ2 or SHORE flat file format,

三、比对命令

这里要分两步走，首先是构建参考基因组的索引，然后才是比对

/home/jmzeng/bio-soft/genomemapper-0.4.4/gmindex \

-i BRCA1.fa -x BRCA1.idx -t BRCA1.meta

首先构建索引，种子长度就用默认的12即可，然后构建完索引如下。

image002

然后进行比对即可

/home/jmzeng/bio-soft/genomemapper-0.4.4/genomemapper \

-i BRCA1.fa -q SRR258835.fastq -M 4 -G 2 -E 4 -o mapped_reads.fl -u unmapped_reads.fl

成功比对的都输出到了mapped_reads.fl -这个文件，未比对上的在unmapped_reads.fl

我有12344条序列，成功比对的只有5276条，但是如果我用精确比对的算法，只有一千五百条是可以比对的，所以用这个允许4个mismatch和2个gap的比对算法，大大提高了比对率。

然后我修改了比对参数可以达到5605,5654,5696的提升。但是没有质的飞跃，估计本身我的这种helicos测序数据错误率就太可怕了。

四，输出结果解读

image004

这个是很规则的tab键分割的文本字符，我就不解读了，大家看readme

六 08

探究各个步骤对snp-calling的影响

Posted on by

做snp-calling时很多标准流程都会提到去除PCR重复这个步骤，但是这个步骤对找snp的影响到底有多大呢？这里我们来探究一下

 

去除PCR重复前 样本名 去除PCR重复后
   106082 BC1-1.snp 103829
   101443 BC1-2.snp 99500
   103937 BC2-1.snp 101833
   102979 BC2-2.snp 101022
   105876 BC3-1.snp 103562
   109168 BC3-2.snp 107052
   107155 BC4-1.snp 104894
   108335 BC4-2.snp 106031
   100236 BC5-1.snp 98417
   102322 BC5-2.snp 100395
   103466 BC6-1.snp 101405
   112940 BC6-2.snp 110611
   113166 BC7-1.snp 110948
   114038 BC7-2.snp 116090
   123670 PC1-1.snp 121697
   111402 PC1-2.snp 109389
   106917 PC2-1.snp 105149
   108724 PC2-2.snp 106776

 

可以看到去除pcr重复这个脚本对snp-calling的结果影响甚小，就是少了那么一千多个snp，脚本如下，我是用picard-tools进行的去除PCR重复，当然也可以用samtools来进行同样的步骤

[shell]

<b>for i in *.sorted.bam</b>

<b>do</b>

<b>echo $i</b>

<b>java  -Xmx120g  -jar /home/jmzeng/snp-calling/resources/apps/picard-tools-1.119/MarkDuplicates.jar \</b>

<b>CREATE_INDEX=true REMOVE_DUPLICATES=True \</b>

<b>ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT METRICS_FILE=/dev/null \</b>

<b>INPUT=$i OUTPUT=${i%%.*}.sort.dedup.bam</b>

<b>done</b>

[/shell]

然后我们首先看看没有产生变化的那些snp信息的改变

head -50  ../rmdup/out/snp/BC1-1.snp  |tail |cut -f 1,2,8

chr1 17222 ADP=428;WT=0;HET=1;HOM=0;NC=0

chr1 17999 ADP=185;WT=0;HET=1;HOM=0;NC=0

chr1 18091 ADP=147;WT=0;HET=1;HOM=0;NC=0

chr1 18200 ADP=278;WT=0;HET=1;HOM=0;NC=0

chr1 24786 ADP=238;WT=0;HET=1;HOM=0;NC=0

chr1 25072 ADP=24;WT=0;HET=1;HOM=0;NC=0

chr1 29256 ADP=44;WT=0;HET=1;HOM=0;NC=0

chr1 29265 ADP=44;WT=0;HET=1;HOM=0;NC=0

chr1 29790 ADP=351;WT=0;HET=1;HOM=0;NC=0

chr1 29939 ADP=109;WT=0;HET=1;HOM=0;NC=0

head -50   BC1-1.snp  |tail |cut -f 1,2,8

chr1 17222 ADP=457;WT=0;HET=1;HOM=0;NC=0

chr1 17999 ADP=196;WT=0;HET=1;HOM=0;NC=0

chr1 18091 ADP=155;WT=0;HET=1;HOM=0;NC=0

chr1 18200 ADP=313;WT=0;HET=1;HOM=0;NC=0

chr1 24786 ADP=254;WT=0;HET=1;HOM=0;NC=0

chr1 25072 ADP=25;WT=0;HET=1;HOM=0;NC=0

chr1 29256 ADP=46;WT=0;HET=1;HOM=0;NC=0

chr1 29265 ADP=46;WT=0;HET=1;HOM=0;NC=0

chr1 29790 ADP=440;WT=0;HET=1;HOM=0;NC=0

chr1 29939 ADP=123;WT=0;HET=1;HOM=0;NC=0

可以看到，同一位点的snp仍然可以找到，仅仅是对测序深度产生了影响

 
然后我们再看看去除PCR重复这个步骤减少了的snp，在原snp里面是怎么样的

perl -alne '{$file++ if eof(ARGV);unless ($file){$hash{"$F[0]_$F[1]"}=1} else {print if not exists $hash{"$F[0]_$F[1]"} } }' ../rmdup/out/snp/BC1-1.snp BC1-1.snp |less

这个脚本就可以把去除PCR重复找到的snp位点在没有去除PCR重复的找到的snp文件里面过滤掉，查看那些去除PCR重复之前独有的snp

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

8.00    8.00   11.00   44.26   25.00 7966.00

图片1

 

可以看到被过滤的snp大多都是测序深度太低了的，如下面的例子

chr1 726325 a 9 CCC.ccc,^:, IEHGHHG/9

chr1 726325 a 5 C.c,^:, IGH/9

 

chr1 726338 g 16 TTT.ttt,,....,,, IHGI:9<HIIFIHC5H

chr1 726338 g 10 T.t,,...,, II:HIIFH5H

 

可以看到这一步还是很有用的，但是怎么说呢，因为最后对snp的过滤本来就包含了一个步骤是对snp的测序深度小于20的给过滤掉

 

但是也有个别的测序深度非常高的snp居然也是被去除PCR重复这个步骤给搞没了！很奇怪，我还在探索之中.

grep 13777 BC1-1.mpileup  |head

chr1 13777 G 263 ........,.C,,,,,.,,,.......,,,..,....,,......,.....c,........,,,,,,,..,...,,,,,.........,......C.......,,,,,,,,,,.....,,,,,,,.,,,..C,,,,,,CC,c,,,...C..,,,,cC.C..CC.CC,,cc,.C...C,,,,CCc,c,,,,,,,c,C.C.CC...C.cc,c...,C.CCcc...,CCC.C.CC..CCC..CC.c,cc,cc,,cc,C.,,^!.^6.^6.^6.^!, HIHIIIIEIEIHGIIIFIHIG?IIIIHIIHIFHIIHICIIIHIIGIEIIGIIIGHIIIIIIHIIHIHIIIIIIIHII1I?GHHHEHHIIEIEHIIEIIHHIIFIIIFHIHIIIIHIHIIHIIHHIIEIIIIIIHIIIIIIIIIG1HIIIIHIHIEHIHIHIIIIIIIIIIIHICIHIIIIIEIIIIHICIHGGIIIIIIIIEHIHIIIIIIHFIGGIIIIGIIIGICIIIHIIIIIIIIIIIHHHIIIIIHIIHDDII>>>>>

grep 13777 BC1-1.rmdup.mpileup  |head

chr1 13777 G 240 ........,.C,,,,,.,,,.......,,,..,....,,......,....c,......,,,,,,,..,...,,,,,.........,......C......,,,,,,,,,,.....,,,,,,,.,,,..C,,,,,,CC,c,,,...C..,,,,cC..CC.CC,cc,.C...C,,,,CCc,c,,,,,,,cC.C.C..C.c,c...,C.CCcc...,CC.C.CCC..C.c,cc,,c,.,,^!.^6.^6.^!, HIHIIIIEIEIHGIIIFIHIG?IIIIHIIHIFHIIHICIIIHIIGIEIIIIIHIIIIIHIIHIHIIIIIIIHII1I?GHHHEHHIIEIEHIIEIHHIIFIIIFHIHIIIIHIHIIHIIHHIIEIIIIIIHIIIIIIIIIG1HIIIIHIHIEHIHIIIIIIIIIIHICIHIIIIIEIIIIHICIHGGIIIIIIHIHIIIIIHFIGGIIIIGIIIGCIIIIIIIIIIHHIIIHIHDII>>>>

 

然后我再搜索了一些

chr8 43092928 . A T . PASS ADP=7966;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:7967:7966:6261:1663:20.9%:0E0:39:39:3647:2614:1224:439

chr8 43092908 . T C . PASS ADP=6968;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:7002:6968:5315:1537:22.06%:0E0:37:38:3022:2293:890:647

chr8 43092898 . T G . PASS ADP=6517;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:6517:6517:4580:1587:24.35%:0E0:38:38:2533:2047:920:667

chr7 100642950 . T C . PASS ADP=770;WT=0;HET=1;HOM=0;NC=0 GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR 0/1:255:771:770:615:155:20.13%:3.9035E-51:38:38:277:338:65:90

终于发现规律啦！！！原来它们的突变率都略高于20%，在没有去处PCR重复之前，是高于snp的阈值的，但是去除PCR重复对该位点的突变率产生了影响，使之未能通过筛选。

 

六 01

Samtools无法同时得到mpileup格式的数据和bcftools格式的数据

Posted on by

 来自于： https://www.biostars.org/p/63429/

I'm using samtools mpileup and would like to generate both a pileup file and a vcf file as output. I can see how to generate one or the other, but not both (unless I run mpileup twice). I suspect I am missing something simple.

Specifically, calling mpileup with the -g or -u flag causes it to compute genotype likelihoods and output a bcf. Leaving these flags off just gives a pileup. Is there any way to get both, without redoing the work of producing the pileup file? Can I get samtools to generate the bcf _from_ the pileup file in some way? Generating the bcf from the bam file, when I already have the pileup, seems wasteful.

Thanks for any help!

我写了脚本来运行，才发现我居然需要两个重复的步骤来得到mpileup格式的数据和bcftools格式的数据，而这很明显的重复并且浪费时间的工作

for i in *sam

do

echo $i

samtools view -bS $i >${i%.*}.bam

samtools sort ${i%.*}.bam ${i%.*}.sorted

samtools index ${i%.*}.sorted.bam

samtools mpileup -f /home/jmzeng/ref-database/hg19.fa  ${i%.*}.sorted.bam  >${i%.*}.mpileup

samtools mpileup -guSDf  /home/jmzeng/ref-database/hg19.fa  ${i%.*}.sorted.bam  | bcftools view -cvNg - > ${i%.*}.vcf

Done

我想得到mpileup格式，是因为后续的varscan等软件需要这个文件来call snp

而得到bcftools格式可以直接用bcftools进行snp-calling

samtools mpileup 命令只有用了-g或者-u那么就只会输出bcf文件

如果想得到mpileup格式的数据，就只能用-f参数。

  • bcftools doesn't work on pileup format data. It works on bcf/vcf files.
  • samtools provides a script called sam2vcf.pl, which works on the output of "samtools pileup". However, this command is deserted in newer versions. The output of "samtools mpileup" does not satisfy the requirement of sam2vcf.pl. You can check the required pileup format on lines 95-99, which is different from output of "samtools mpileup".

 

五 29

用R语言的RCurl包结合XML包批量下载生信课件

Posted on by

首先是宾夕法尼亚州立大学（The Pennsylvania State University缩写PSU)的生信课件下载，这个生信不仅有课件，而且在中国的优酷视频网站里面还有全套授课视频，非常棒！

image001

课程主页是http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html

可以看出所有的课件pdf链接都在这一个页面，所以是非常简单的代码！

下面是R代码：

library(XML)

library(RCurl)

library(dplyr)

psu_edu_url='http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html';

wp=getURL(psu_edu_url)

base='http://www.personal.psu.edu/iua1/courses/file';

#pse_edu_links=getHTMLLinks(psu_edu_url)

psu_edu_links=getHTMLLinks(wp)

psu_edu_pdf=psu_edu_links[grepl(".pdf$",psu_edu_links,perl=T)]

for (pdf in psu_edu_pdf){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

因为这三十个课件都是接近于10M，所以下载还是蛮耗时间的

image003

其实R语言里面有这个down_url函数，可以直接下载download.file(down_url,filename)

然后我开始下载德国自由大学的生信课件，这次不同于宾夕法尼亚州立大学的区别是，课程主页里面是各个课题的链接，而pdf讲义在各个课题里面，所以我把pdf下载写成了一个函数对我们的课题进行批量处理

library(XML)

library(RCurl)

library(dplyr)

base="http://www.mi.fu-berlin.de/w/ABI/Genomics12";

down_pdf=function(url){

links=getHTMLLinks(url)

pdf_links=links[grepl(".pdf$",links,perl=T)]

for (pdf in pdf_links){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

}

down_pdf(base)

list_lecture= paste("http://www.mi.fu-berlin.de/w/ABI/GenomicsLecture",1:15,"Materials",sep="")

for ( url in list_lecture ){

cat("Now we process the ",url ,"\n")

try(down_pdf(url))

}

image005

同样也是很多pdf需要下载

接下来下载Minnesota大学的关于生物信息的教程的ppt合集

主页是： https://www.msi.umn.edu/tutorial-materials

 

这个网页里面有64篇pdf格式的ppt，还有几个压缩包，本来是准备写爬虫来爬去的，但是后来想了想有点麻烦，而且还不一定会看，反正也是玩玩

就用linux的命令行简单实现了这个爬虫功能。

curl https://www.msi.umn.edu/tutorial-materials >tmp.txt

perl -alne '{/(https.*?pdf)/;print $1 if $1}' tmp.txt >pdf.address

perl -alne '{/(https.*?txt)/;print $1 if $1}' tmp.txt

perl -alne '{/(https.*?zip)/;print $1 if $1}' tmp.txt >zip.address

wget -i pdf.address

wget -i pdf.zip

这样就可以啦！

 

用爬虫也就是几句话的事情，因为我已经写好啦下载函数，只需要换一个主页即可下载页面所有的pdf文件啦！

 

五 29

R语言对Ozone数据处理笔记

Posted on by

一．首先加载一些包，这样才能获得该ozone数据

数据集介绍：

Ozone数据集是一个三维数组，记录了24×24个空间网格内，从1995年1月到2000年12月，共72个时间点上，中美洲每月的平均臭氧水平。

前两维分别表示纬度和经度，第三维表示时间。

加载包的代码如下：

library("MASS", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("ggplot2", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("plyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("dplyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("reshape2", lib.loc="C:/Program Files/R/R-3.1.1/library")

二．我们首先简单看看第一个地点的72个月的臭氧水平变化图。

plot(1:72,ozone[1,1,],type="l")

box(lty = '1373', col = 'red')

grid(nx=NA,ny=NULL,lty=1,lwd=1,col="gray")

看起来还算是蛮有规律的。

image001

三．然后把这72个月的数据分成年份来画图

绘图第一种方式如下：

value <-ozone[1, 1, ]

plot(value[1:12],type="b",pch=19,lwd=2,xaxt="n",col="black",

xlab="month",ylab="value")

axis(1,at=1:12,labels=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))

lines(value[13:24],col="red",type="b",pch=19,lwd=2)

lines(value[25:36],col="orange",type="b",pch=19,lwd=2)

lines(value[37:48],col="purple",type="b",pch=19,lwd=2)

lines(value[49:60],col="blue",type="b",pch=19,lwd=2)

lines(value[61:72],col="green",type="b",pch=19,lwd=2)

legend("bottomright",legend=c("1995","1996","1997","1998","1999","2000"),

lty=1,lwd=2,pch=rep(19,6),col=c("black","red","orange","purple","blue","green"),

ncol=1,bty="n",cex=1.2,

text.col=c("black","red","orange","purple","blue","green"),inset=0.01)

是首先画第一年的，然后逐年添加一条线，然后画图例，算是蛮复杂的。

image003

还有一个简单的方法，就是用ggplot这个包来画。

values <-ozone[1, 1, ]

months=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

years=c("1995","1996","1997","1998","1999","2000")

dat=data.frame(month=rep(months,6),value=values,year=rep(years,each=12))

ggplot(dat,aes(x=month,y=value,group=year,colour=year))+geom_line()

image005

四：测试一下稳健回归模型

稳健回归是加权最小二乘回归，或称文艺最小二乘回归。

MASS 包中的 rlm命令提供了不同形式的稳健回归拟合方式。

回归分析就是用数理统计的方法，研究自然界中变量之间存在的非确定的相互依赖和制约关系，并把这种非确定的相互依赖和制约关系用数学表达式表达出来。其目的在于利用这些数学表达式以及对这些表达式的精度估计，对未知变量作出预测或检验其变化，为决策服务。

介绍几个线性回归（linear regression）中的术语：

残差（Residual）: 基于回归方程的预测值与观测值的差。

离群点（Outlier）: 线性回归（linear regression）中的离群点是指对应残差较大的观测值。也就是说，当某个观测值与基于回归方程的预测值相差较大时，该观测值即可视为离群点。 离群点的出现一般是因为样本自身较为特殊或者数据录入错误导致的，当然也可能是其他问题。

杠杆率（Leverage）: 当某个观测值所对应的预测值为极端值时，该观测值称为高杠杆率点。杠杆率衡量的是独立变量对自身均值的偏异程度。高杠杆率的观测值对于回归方程的参数有重大影响。

影响力点：（Influence）: 若某观测值的剔除与否，对回归方程的系数估计有显著相应，则该观测值是具有影响力的，称为影响力点。影响力是高杠杆率和离群情况引起的。

Cook距离（Cook's distance）: 综合了杠杆率信息和残差信息的统计量。

values <-ozone[1, 1, ]

month.abbr=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

month_72 <-factor(rep(month.abbr, length = 72), levels = month.abbr)

deseas1 <-rlm(value ~ month_72 -1)

summary(deseas1)

image007

五．对我们的24*24个地理位置的数据都做稳健回归分析

deseasf<-function(value) rlm(value ~ month_72 -1, maxit= 50)

models <-alply(ozone, 1:2, deseasf) #model是一个list，储存着24*24次稳健回归的结果

failed <-laply(models, function(x) !x$converged)

coefs<-laply(models, coef)#coefs是一个三维数组，记录了所有24×24个位置中每个位置的12个系数

dimnames(coefs)[[3]] <-month.abbr

names(dimnames(coefs))[3] <-"month"

deseas<-laply(models, resid) #deseas是一个三维数组，记录了所有24×24个位置中每个位置的72个残差

dimnames(deseas)[[3]] <-1:72

names(dimnames(deseas))[3] <-"time"

六．对我们的稳健回归系数的三维矩阵进行降维处理，方便画图

通过reshape包可以对三维数组进行降维

coefs_df<-melt(coefs)

head(coefs_df)

lat   long month   value

1 -21.2 -113.8   Jan 264.3964

2 -18.7 -113.8   Jan 261.3284

3 -16.2 -113.8   Jan 260.9643

4 -13.7 -113.8   Jan 258.9999

5 -11.2 -113.8   Jan 255.9999

6 -8.7 -113.8   Jan 254.9999

可以看到第三维的month成功被降维了

还可以通过plyr这个数据工厂包来进行降维

coefs_df<-ddply(coefs_df, .(lat, long), transform, avg= mean(value), std= value / max(value))

>head(coefs_df)   lat   long month   value     avg       std1 -21.2 -113.8   Jan 264.3964 268.6604 0.92770682 -21.2 -113.8   Feb 259.2036 268.6604 0.90948623 -21.2 -113.8   Mar 255.0000 268.6604 0.89473684 -21.2 -113.8   Apr 252.0052 268.6604 0.88422885 -21.2 -113.8   May 258.5089 268.6604 0.90704866 -21.2 -113.8   Jun 265.3387 268.6604 0.9310129

可以看到，不仅成功降维了，还添加了几个属性变量

 

七．最后对降维的coef系数数据画热图

 

coef_limits<-range(coefs_df$value)

coef_mid<-mean(coefs_df$value)

monthsurface<-function(mon)

{

df<-subset(coefs_df, month == mon)

qplot(long, lat, data = df, fill = value, geom="tile") +

scale_fill_gradient(limits = coef_limits,

low = "lightskyblue", high = "yellow")

}

pdf("ozone-animation.pdf", width = 8, height = 8)

l_ply(month.abbr, monthsurface, .print = TRUE)

dev.off()

会在当前R的工作目录下面看到一个pdf文件，里面储存着12个月的在24*24个地理位置的系数热图。

image009

五 28

org.Xx.eg.db系列包概述

Posted on by

在bioconductor的官网里面可以查到共有111个系列包，基本上跨越了我们常见的物种啦！

org.Xx.eg.db系列包介绍65

斑马鱼：Bioconductor - org.Dr.eg.db - /packages/release/data/annotation/html/org.Dr.eg.db.html

Details biocViews AnnotationData , Danio_rerio , OrgDb Version 3

拟南芥：Bioconductor - org.At.tair.db - /packages/release/data/annotation/html/org.At.tair.db.html

Details biocViews AnnotationData , Arabidopsis_thaliana , OrgDb Version 3

小鼠：Bioconductor - org.Mm.eg.db - /packages/release/data/annotation/html/org.Mm.eg.db.html

Details biocViews AnnotationData , Mus_musculus , OrgDb , mouseLLMappings Version 3

人类：Bioconductor - org.Hs.eg.db - /packages/release/data/annotation/html/org.Hs.eg.db.html

Details biocViews AnnotationData , Homo_sapiens , OrgDb , humanLLMappings Version 3

对这些系列包的函数都一样，包括以下几个：

columns(x)  keytypes(x)  keys(x, keytype, ...)  select(x, keys, columns, keytype, ...)  saveDb(x, file)  loadDb(file, dbType, dbPackage, ...)

 

 

这些包就是bioconductor已经做好的数据库，我们可以根据定义好的ID号来进行任意的基因转换，现在支持的信息有一下几种！

keytypes(org.Hs.eg.db)

[1] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"      "ACCNUM"       "ALIAS"        "CHR"

[8] "CHRLOC"       "CHRLOCEND"    "ENZYME"       "MAP"          "PATH"         "PMID"         "REFSEQ"

[15] "SYMBOL"       "UNIGENE"      "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "GENENAME"     "UNIPROT"

[22] "GO"           "EVIDENCE"     "ONTOLOGY"     "GOALL"        "EVIDENCEALL"  "ONTOLOGYALL"  "OMIM"

[29] "UCSCKG"

这些包的确非常有用，大家可以看我博客里面关于它们的介绍！！！

 

五 28

菜鸟团第二次作业的部分答案

Posted on by

> library(org.Hs.eg.db)

载入需要的程辑包：AnnotationDbi载入需要的程辑包：stats4载入需要的程辑包：GenomeInfoDb载入需要的程辑包：S4Vectors载入需要的程辑包：IRanges载入程辑包：‘AnnotationDbi’The following object is masked from ‘package:GenomeInfoDb’:     species载入需要的程辑包：DBI

 

1、人共有多少个entrez id的基因呢？

x <- org.Hs.egENSEMBLTRANS

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(x)

[1] 47721

可知共有47721个基因都是有entrez ID号的

2、能对应转录本ID的基因有多少个呢？

length(xx)

[1] 20592

可以看到共有20592个基因都是有转录本的！

2、能对应ensembl的gene ID的基因有多少个呢？

x <- org.Hs.egENSEMBL

# Get the entrez gene IDs that are mapped to an Ensembl ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 26019

可以看到只有26019是有ensembl的gene ID的

3、那么基因对应的转录本分布情况如何呢？

table(unlist(lapply(xx,length)))

菜鸟团第二次作业的部分答案863

可以看出绝大部分的基因都是20个转录本一下的，但也有极个别基因居然有高达两百个转录本，很可怕！

4、那么基因在染色体的分布情况如何呢？

x <- org.Hs.egCHR

# Get the entrez gene identifiers that are mapped to a chromosome

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

> length(x)

[1] 47721

> length(xx)

[1] 47232

可以看到有接近五百个基因居然是没有染色体定位信息的！！！

table(unlist(xx))

用barplot函数可视化一下，如图

 

菜鸟团第二次作业的部分答案1209

6、那么有多多少基因是有GO注释的呢？

x <- org.Hs.egGO

# Get the entrez gene identifiers that are mapped to a GO ID

mapped_genes <- mappedkeys(x)

# Convert to a list

xx <- as.list(x[mapped_genes])

length(xx)

[1] 18229

> length(x)

[1] 47721

可以看到只有18229个基因是有go注释信息的。

那么基因被注释的go的分布如何呢？

菜鸟团第二次作业的部分答案1477

可以看到大部分的基因都是只有30个go的，但是某些基因特别活跃，高达197个go注释。

还有kegg和omin数据库的我就不写了！

五 28

实战R语言bioconductor的seqinr包探究人的所有转录本的性质

Posted on by

首先安装这个包

source("http://bioconductor.org/biocLite.R")

biocLite("seqinr")

然后加载包，并读取我们的CDS.fa文件

library("seqinr")

human_cds=read.fasta("CDS.fa")

#这一个步骤非常耗时间，可能是因为我们的转录本文件有十万多个转录本的原因吧

str(human_cds) #查看可知读入了一个list，其中每个转录本都是list的一个元素

List of 100778

$ ENST00000415118:Class 'SeqFastadna'  atomic [1:8] g a a a ...

.. ..- attr(*, "name")= chr "ENST00000415118"

.. ..- attr(*, "Annot")= chr ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene tran"| __truncated__

$ ENST00000448914:Class 'SeqFastadna'  atomic [1:13] a c t g ...

.. ..- attr(*, "name")= chr "ENST00000448914"

.. ..- attr(*, "Annot")= chr ">ENST00000448914 havana_ig_gene:known chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene tran"| __truncated__

对list的每个元素都有几种函数可以处理得到信息：

Length,table,GC,count

其中count函数很有趣，数一数序列里面的这些组合出现的次数

count(dengueseq, 1)

count(dengueseq, 2)接下来我们随机取human_cds这个list的一个元素用这几个函数对它处理一下

> tmp=human_cds[[1]]

> tmp

[1] "g" "a" "a" "a" "t" "a" "g" "t"

attr(,"name")

[1] "ENST00000415118"

attr(,"Annot")

[1] ">ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene"

attr(,"class")

[1] "SeqFastadna"

再看看函数的结果

> length(tmp)

[1] 8

> table(tmp)

tmp

a g t

4 2 2

> GC(tmp)

[1] 0.25

> count(tmp,1)

 

a c g t

4 0 2 2

> count(tmp,2)

 

aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt

2  0  1  1  0  0  0  0  1  0  0  1  1  0  0  0

>

还是挺好用的，接下来我们应用R的知识来对着十万多个转录本进行一些简单的总结

human_cds_length=unlist(lapply(human_cds,length))

human_cds_gc=unlist(lapply(human_cds,GC))

这样就得到了所有转录本的长度和GC含量信息

然后我们简单统计一下，并画几个图表吧！

> summary(human_cds_length)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

3     366     699    1132    1425  108000

> summary(human_cds_gc)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

0.1467  0.4577  0.5285  0.5264  0.5932  0.8917

可以看到还是有很多很极端的转录本的存在！

最长的转录本也不过10k，而我记得最长的基因高达8M，看了内含子远大于外显子呀。

但是GC含量有很多高于80%，这些基因在二代测序的研究中是一个盲区。

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2075

这些极端基因可以通过biomaRt等包进行注释，得到gene名和功能信息。

 

hist(human_cds_gc)

hist(log10(human_cds_length))

GC含量分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2177

长度分布如图

实战R语言bioconductor的seqinr包探究人的所有转录本的性质2186

附表：

http://www.bioinformatics.org/sms/iupac.html 所有字符的碱基氨基酸意义表格

 

五 25

用R画GO注释二级分类统计图

Posted on by

群里有朋友问这个图怎么画，我想了想，这肯定是ggplot完成的，非常简单，但是菜鸟们缺乏实践，可能会困惑，所以我模拟数据画了一个！

图片1

首先构造数据

dat=data.frame(name=LETTERS[1:21],
 number=abs(rnorm(21)*10),
 type=c(rep("BP",7),rep("CC",7),rep("MF",7))
)
# 请务必自己查看dat是一个什么数据，print出来即可
# 然后对这个数据画图，一行代码即可！！！
library(ggplot2)
ggplot(dat,aes(x=name,y=number,fill=type))+geom_bar(stat="identity")+coord_flip()

看起来是不是很像回事啦！细节我就懒得调控啦！

图片2

其实自己搜索即可！坐标轴和主题都是可以控制的

http://rstudio-pubs-static.s3.amazonaws.com/3364_d1a578f521174152b46b19d0c83cbe7e.html

http://docs.ggplot2.org/0.9.3.1/coord_flip.html

​
五 21

Biostrings包简介

Posted on by

首先讲讲它的对象

有下面几个字符串对象BString, DNAString, RNAString and AAString可以通过以下代码构造它们：

b <- BString("I am a BString object")

d <- DNAString("TTGAAAA-CTC-N")

这两个对象的区别是DNAstring对象对字符串的要求严格一些，只有IUPAC字符和+-字符可以。

对构造好的对象可以通过下标来取子字符串对象，也可以通过subseq来取，但是子字符串仍然是数据对象，只有通过toString函数才能把它们转化成字符串。

用length(dd2)和nchar(toString(dd2))都可以找到我们Biostrings对象的长度。但是后者速度会很慢。

Views(RNAString("AU"), start=0, end=2)这个函数可以把string对象任意截取成list

start, end and width可以作用于我们截取的list，判断list里面的元素在原来的string对象上面的起始终止及长度信息。

 

接下来讲这个包带有的一个比对函数！

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] succ--eed subject: [1] supersede score: -33.99738

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",type = "local")

Local PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su subject: [1] su score: 5.578203

> pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",gapOpening = 0, gapExtension = 1)

Global PairwiseAlignmentsSingleSubject (1 of 2)pattern: [1] su-cce--ed subject: [1] sup--

可以看出这个比对函数可以调整的参数实在是太多了，而且改变参数之后比对情况大不一样，还有很多参数就不一一细讲了。

这个比对结果可以赋值给一个变量，保存比对的对象。

psa1 <- pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede")

class(psa1)

summary(psa1)

class(pattern(psa1))

class(summary(psa1))

score(psa2)

还可以自己构建打分矩阵来进行比对。

submat <-

+ matrix(-1, nrow = 26, ncol = 26, dimnames = list(letters, letters))

diag(submat) <- 0

Biostrings包简介1454

psa2 <-pairwiseAlignment(pattern = c("succeed", "precede"), subject = "supersede",substitutionMatrix = submat, gapOpening = 0, gapExtension = 1)

我们的包还自带了两个非常流行的氨基酸比对矩阵PAM和BLOSUM

ls("package:Biostrings")可以查看这个包所有的对象。

data(package="Biostrings")可以查看这个包所有的数据对象

还有很多其它函数

还可以去除adaptor，挺好玩的

既然有配对比对函数，那么就有多重比对函数！

我们可以读取clustaW, Phylip and Stolkholm这几种不同的比对结果文件来构造多重比对对象。

library(Biostrings)这个包里面自带了两个文件，我们可以示范一下构建对象。

origMAlign <- readDNAMultipleAlignment(filepath = system.file("extdata", "msx2_mRNA.aln", package="Biostrings"), format="clustal")

phylipMAlign <- readAAMultipleAlignment(filepath = system.file("extdata","Phylip.txt", package="Biostrings"),format="phylip")

 

对构造好的多重比对对象就可以构建进化树啦，代码如下！

sdist <- stringDist(as(origMAlign,"DNAStringSet"), method="hamming")

> clust <- hclust(sdist, method = "single")

> pdf(file="badTree.pdf")

> plot(clust)

> dev.off()

Biostrings包简介2345

五 21

Bioconductor的DO.db包介绍

Posted on by

Bioconductor的包都是同样的安装方法：

source("http://bioconductor.org/biocLite.R");biocLite("DO.db")

还有GO.bd包是完全一模一样的规则！！！

加载这个包可以发现它依赖于好几个其它的包，这也是我比较喜欢R的原因，它会自动把它需要的包全部安装加载进来，不需要自己一个个调试！

> library(DO.db)

载入需要的程辑包：AnnotationDbi

载入需要的程辑包：stats4

载入需要的程辑包：GenomeInfoDb

载入需要的程辑包：S4Vectors

载入需要的程辑包：IRanges

载入需要的程辑包：DBI

> help(DO.db)

> ls("package:DO.db")

[1] "DO"          "DO_dbconn"   "DO_dbfile"   "DO_dbInfo"   "DO_dbschema" "DOANCESTOR"  "DOCHILDREN"  "DOID"        "DOMAPCOUNTS"

[10] "DOOBSOLETE"  "DOOFFSPRING" "DOPARENTS"   "DOSYNONYM"   "DOTERM"      "DOTerms"     "Secondary"   "show"        "Synonym"

[19] "Term"

这个包里面有19个数据对象！都是比较高级的S4对象。

比如我们可以拿DOTERM[1:10]这个小的数据对象来做例子！example=DOTERM[1:10]

因为example是一个高级对象，所以无法直接查看，需要用as.list方法来查看

> as.list(example)

$`DOID:0001816`DOID: DOID:0001816Term: angiosarcomaSynonym: DOID:267Synonym: DOID:4508Synonym: "hemangiosarcoma" EXACT []Secondary: DOID:267Secondary: DOID:4508

~~~~~~~~~~~~共十个DO条目

对每一个DO条目来说都有DOID,Term,Synony这些函数可以取对应的值。

下面是对DO的有向无环图的数据解读

xx <- as.list(DOANCESTOR)可以查看每个DO与它所对应的上级条目DO，每个DO都会有不止一个的上级DO。

xx <- as.list(DOPARENTS)可以查看每个DO与它所对应的父条目DO，每个DO都有且只有一个父DO。

xx <- as.list(DOOFFSPRING)可以查看每个DO与它所对应的下级DO的关系列表，大多数DO都不止一个子条目DO，所有的下级DO都会列出。

xx <- as.list(DOCHILDREN)以查看每个DO与它所对应的子条目DO的关系列表，大多数DO都不止一个子条目DO。

还有Lkeys(DOTERM)可以查看数据库里面的所有的DO条目的ID号

> head(keys(DOTERM))

[1] "DOID:0000000" "DOID:0001816" "DOID:0002116" "DOID:0014667" "DOID:0050004" "DOID:0050012"

dbmeta(GO_dbconn(), "GOSOURCEDATE")

可以查看这个DO库的制备时间

> dbmeta(DO_dbconn(), "DOSOURCEDATE")

[1] "20140417"