Monthly Archives: 7月 2015

七 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