科学网—标签

相关帖子	版块	作者	回复/查看	最后发表

mashengwei 2018-1-25 22:39

call variants from wheat RNA_seq 1 24 本期作者：Neal call variants from wheat RNA_seq 上周我们推送的何中虎研究员的文章（中国小麦产业发展与科技进步——小麦里我见过的最豪华作者阵容），目前阅读人数是8400多人，也是迄今为止我们公众号阅读人数最高的推送之一，这也直接让我们的关注人数猛增到3000人。过去一段时间，有不少小麦育种老师和专家也关注了我们，我们非常欢迎各位分享育种方面的经验。今天要说一些说从RNA seq数据里得到序列变异的步骤。首先要交代一下背景，我们要研究一个EMS突变体（已回交多次），前期已经将突变基因定位到一个染色体区间，根据中国春参考序列，我们已经发现一个候选基因已经在水稻里被报道过了，测序发现，该基因确实发生了变异。后期安排了RNA seq实验，想在机制方面做一个有益的尝试。另外我们也想通过这样一个RNA_seq数据，比较在定位区间内还有那些基因发生了变异，因此就有了今天这样一个推送。这里还要特别强调一点，这不是混池数据，分析过程中的一些参数请根据实验目的调整。使用STAR将reads mapping至小麦基因组，然后使用sentieon流程（本质是GATK）call variant，接着使用SnpSift筛选高质量SNP，结合EMS诱变的特点，进一步排除可能的假阳性SNP，最后获得大概300个SNP，使用SnpEff注释SNP。根据前期遗传定位结果，我们发现只有4个SNP位于我们的区间内(20Mb)，但是只有一个SNP导致蛋白提前终止，该SNP所在的基因其在水稻里的直系同源基因已被报道，突变之后与我们的突变体表型非常相似。下面是具体的流程，如果有兴趣欢迎交流。这个需要具有一定的高通量数据分析基础，其他的就没有什么特别的地方了。 # 工作目录 cd / data / rna_seq / genome / #构建 A 基因组 index STAR -- runThreadN 10 -- runMode genomeGenerate -- genomeDir ./ -- genomeFastaFiles CS_A_genome_part . fasta -- sjdbFileChrStartEnd TGACv1_part_A . ss -- limitGenomeGenerateRAM 68800833920 #!/usr/bin/env python # -*- coding: utf-8 -*- __author__ = 'wheatomics' import subprocess # inpu.txt放着fastq文件的名字 with open ( 'input.txt' , 'r' ) as f : for line in f : line = line . strip (). split () fq1 , fq2 = line print fq1 , fq2 # 1. Mapping reads with STAR proc = subprocess . Popen ( , 'SM:' + fq1 . split ( '_' ) , 'PL:ILLUMINA' ], shell = False ) proc . wait () proc = subprocess . Popen ( , shell = False ) proc . wait () proc = subprocess . Popen ( , shell = False ) proc . wait () proc = subprocess . Popen ( , shell = False ) proc . wait () # 2. Metrics proc = subprocess . Popen ( , shell = False ) proc . wait () proc = subprocess . Popen ( + '-metrics-report.pdf' , 'gc=gc_metrics.txt' , 'qd=qd_metrics.txt' , 'mq=mq_metrics.txt' , 'isize=is_metrics.txt' ], shell = False ) proc . wait () # 3. Remove Duplicate Reads proc = subprocess . Popen ( , shell = False ) proc . wait () proc = subprocess . Popen ( , shell = False ) proc . wait () # 4. Split reads at Junction proc = subprocess . Popen ( , shell = False ) proc . wait () # 5. Indel realigner proc = subprocess . Popen ( + '.realigned.bam' ], shell = False ) proc . wait () Call SNP #此处只统计了unique mapped的reads sentieon driver - r / data2 / masw_data / rna_seq / dqyRNA - seq / masw_analysis / genome / CS_A_genome_part . fasta -- read_filter MapQualFilter , min_map_qual = 60 - t 10 - i WT . realigned . bam - i br . realigned . bam -- algo Genotyper -- emit_conf 20 -- call_conf 20 WT_br_UG . vcf 筛选SNP # # 要注意，EMS诱变的变异一般是C/T和A/G的变异，其他类型的变异频率很低很低。另外，突变里的关键功能变异理论上与参考序列不同。EMS mutations result in G:C to A:T mutations, whereas false positives could be any change. Thus, we retained only the alleles that corresponded to G:C to A:T mutations using SnpSift cat WT_ms_UG . vcf | java - jar / data / snpEff / SnpSift . jar filter (QUAL 30) (MQ 40) (QD 5) (FS 30.0) GEN .DP 5 GEN .DP 5 (((GEN .GT = '0/0') (GEN .GT = '1/1')) | ((GEN .GT = '1/1') (GEN .GT = '0/0'))) WT_ms_UG_filtered . vcf #上述筛选参数不是固定的，要根据实验和分析结果调整。具体每个参数表示的意思，请Google搜索SnpSift即可。 vcf文件里的一条染色体还是拆开的，需要合并成一个整体。 #!/usr/bin/env python # -*- coding: utf-8 -*- __author__ = 'wheatomics' chr = , , , , , , , , , , , , , , , , , , , , ] with open ( 'WT_br_UG_filtered.vcf' , 'r' ) as f : for line in f : if line . startswith ( '#' ): print line , else : line = line . replace ( '_part1' , '' ) line = line . strip (). split ( '' ) if line . endswith ( 'part2' ): for i in chr : if line . split ( '_' ) == i : line = int ( line ) + int ( i ) line = line . split ( '_' ) for m in line : print str ( m ) + '' , print line + '' , 最好再修改下vcf的表头信息 ##contig=ID=chr1A,length=594102056,assembly=unknown ##contig=ID=chr2A,length=780798557,assembly=unknown ##contig=ID=chr3A,length=750843639,assembly=unknown ##contig=ID=chr4A,length=744588157,assembly=unknown ##contig=ID=chr5A,length=709773743,assembly=unknown ##contig=ID=chr6A,length=618079260,assembly=unknown ##contig=ID=chr7A,length=736706236,assembly=unknown ##contig=ID=chrUn,length=480980714,assembly=unknown 转换之后就要统计SNP的信息，比如染色体上的SNP个数等 #使用snpEff注释SNP java - Xmx8g - jar snpEff . jar IWGSCv1 . 0 WT_itr_UG_filtered_whole . vcf WT_itr_UG_filtered_whole_eff . vcf 根据上述结果，即可进行下一步的分析。 SnpEff这个我们前面也介绍过，可以参考使用SnpEff 对SNP结果进行分析。欢迎关注 “ 小麦研究联盟 ”，了解小麦新进展投稿、转载、合作以及信息分布等请联系： wheatgenome 原始链接

4353 次阅读|0 个评论

请不要用EMS给我发快件

孙军 2014-4-21 15:31

看看这速度，973的签字快件，所以你就明白为什么这么多人选择私营的快递。我们这个学校在天津偏远的经济技术开发区的最边缘这个现象尤其明显，“天高皇帝远”，这些国营有固定位置的人更慢。经常下午2点半就不能寄送了。我推荐使用顺丰。邮件号码：ER696503595CS 您的邮件于 2014-04-21 09:40:00 （天津市邮政速递物流公司滨海新区分公司黄海路揽投部）投递并签收投递结果：李*代收处理时间处理地点邮件状态 2014-04-18 12:05:00 杭州市文三路数码市场速递经营部收寄 2014-04-18 12:05:00 杭州市文三路数码市场速递经营部揽收 2014-04-18 15:37:44 杭州市文三路数码市场速递经营部离开处理中心,发往杭州市 2014-04-18 17:22:50 杭州市离开处理中心,发往天津市 2014-04-19 10:52:09 天津市离开处理中心,发往天津市邮政速递物流公司滨海新区分公司黄海路揽投部 2014-04-19 15:40:00 天津市邮政速递物流公司滨海新区分公司黄海路揽投部安排投递 2014-04-19 16:00:00 天津市邮政速递物流公司滨海新区分公司黄海路揽投部未妥投 2014-04-19 16:05:16 天津市邮政速递物流公司滨海新区分公司黄海路揽投部到达处理中心 2014-04-21 08:58:00 天津市邮政速递物流公司滨海新区分公司黄海路揽投部安排投递 2014-04-21 09:40:00 天津市邮政速递物流公司滨海新区分公司黄海路揽投部投递并签收

5018 次阅读|0 个评论

EMS，想说爱你不容易

lixuekuan 2010-6-2 09:53

5月28日上午9点46分使用EMS从太原发往陕西延长炼化公司一封公函，到现在已经整整5天（120个小时）了，收件人到现在没有看到信函。从网上查找，没有收到的原因竟然是收件人不在，真会胡说八道，我们电话联系，收件人就在单位急切盼望收到此信函。从离开太原到西安竟然走了50个小时，看来这不是快递，开拖拉机也到了。从洛川县到交口河镇只10公里走了28个小时，昨天下午到达交口河镇至今没有到收件人手里，不知道这还是不是快递。 EMS，官办快递，想说爱你不容易。

个人分类: 生活点滴|4332 次阅读|12 个评论

我和火车不得不说的故事(3)

lxwang 2008-11-22 14:49

那次联系到了火车票，却迟迟拿不到火车票，这就又牵出另外一件事：好不容易买了两张高价票，让上海的朋友给发EMS邮寄过来，本以为EMS有多快捷方便，有多安全可靠，可是我的经历却让我很失望。邮件是前一天下午发出，而第二天一直等到下午也没有收到，又不是很远的地方，从上海到无锡，不是有所谓的次晨达吗！等到第二天下午，我在网上查到，邮件是早上8点发到无锡的，无锡 10点多发往所在地的邮政支局，然后就没有了下文。我打11185，人工服务，一个普通话说得不是很标准的女的大大咧咧的告诉我，10点发往所在地邮政支局，至少下午3，4点才能到，那时候不可能给你送了，只能等到第二天。当时我也没说什么，只是想，这无锡的一个区发过来比上海发过来还慢啊，既然是第二天，那就等吧。第二天，本来熬了一夜做试验的我丝毫没有睡意，等着手机响，等到11点多仍然没有回应。我上网一查询，前一天11点就投妥，我在字面上理解，估计是妥善投递的意思吧，代收人写着一个人的名字。我想可能是所在单位的人帮我收了，所以我去找了所里办公室的人问这个人怎么联系，可人家说全单位就没有这么一个人。我这下慌神了，要知道邮件里面是我两三天回家就要用的高价票啊！我又开始打11185，可打了十几遍都没有打通，老是说业务忙，请稍等，打到最后，听到等候的那个音乐我都头疼了，我才放弃，我怀疑是不是有人在上班啊。我想到打114，查询当地那个邮政支局的联系电话，好在查到了。电话那边一个女的告诉我，邮件是所里专门负责拿信的那个男人的老婆去拿的。听到这句，我又开始联系所里负责拿信的那个人，可我打听到人家老婆也不叫替我代收邮件的那名字啊，可我也没放弃努力，千方百计联系那个人，碰巧办公室都没有那个人的联系方式，办公室没有问到，我就打电话给熟悉的人问也说不知道，当天是周末他也没有来上班，那个下午我都不知道是怎么过来的，焦急等待，无所事事，什么都不想干，满脑子都是要是别人拿了我找不到怎么办啊。下午去食堂吃饭，我问食堂的服务员，我问招待所前台的服务员，后来从大门门卫那得知有那个人的联系方式，让我吃完饭晚上去拿，毕竟人家也要吃饭。我耐住性子吃了一点饭，吃完饭就去了门卫那，可能人家回家吃饭去了，不见人，等了一会还不来，就回宿舍等了一会，又去还是不见人，心里那个急啊，就在往回走的时候才碰到人家，又走到门卫室才拿到那个人的电话。回来打通电话，电话那头一个小女孩声音甜甜的接电话，我就说找你爸爸，一个男人接过了电话，可能是我比较着急，可我知道我还是很诚恳的问人家，那人态度很恶劣，说绝对不可能，他要拿到邮件肯定是当天就会投递。我还是很客气的说我打电话到邮局，说是你夫人去拿的啊，要不你问问，人家还是一口咬定不可能有，就粗暴的挂了电话，让我怀疑刚才那小女孩是不是他女儿啊！我刚挂了电话，就有电话打进，接通，原来是那人的老婆，她问了情况，我说是11点半左右有人代收的，她说她是早上8点去拿的，拿回来三个也没有我的，而11点多签收的就肯定不是她拿回来的，她的解释还说得过去，至少比起刚才那个男的，也就是她丈夫态度好多了。这时候已经是晚上了，我又打11185，终于有人接了，我说我的EMS有人代收了，可是我根本找不到那个人！人家慢条斯理的告诉我明天早上打一个负责投递的电话查询，就这样挂了。一个焦急等待、担心的夜晚，第二天一早，打电话，终于联系到，说我的那个邮件收件人姓名不清楚，电话也不清楚，现在还在邮局。我心里就纳闷了，这样也算投妥？可还是比较欣喜的，毕竟找到了，终于可以回家了！我当即决定，我自己去邮局拿。当时可是飘着大雪啊，南方这么大的雪是很少见的，而气温比较低，雪下到地上，很快在地上就化了，门前正在修路，一有水就是烂泥，怕一时打不到的，我就骑着自行车，冒着大雪往那个邮局赶，雪飘到我的脸上，我的眼睛里，我想可能是逆风，回来时就好了，可是后来回来的时候，情形还是一样，雪直往我脸上吹，我一想，可能是雪还是往下飘的，只是我在往前走，所以无论是去还是回来的时候，我都在往前走，所以就感觉雪是从四面八方吹到我脸上。走了很远才找到那个邮局，拿到邮件一看傻眼了，其实也不怪人家投递人员，在收件人姓名，电话一栏里，什么都没有，还不是看不清楚，是什么都没有！！！这是怎么回事啊，人家给我解释，可能是当时填的时候比较轻，没有印到这一层上。我能怪谁？想想应该是那个收件的邮局工作人的问题，他收件的时候，这样也能行？他为什么不检查一下！这是EMS，号称最快捷，最安全的特快专递，邮局人员如此服务，如此不负责任，怎么行！哎，什么也不说了，至少现在拿到了，至少可以回家了。这就是我的EMS，我的特快专递，我自己冒大雪去邮局拿回来，从上海到无锡已经是第四天才拿到的EMS!!!

个人分类: 人生经历，智慧感悟|2439 次阅读|0 个评论

帐号		自动登录	找回密码
密码			注册

关闭安全验证

标签: EMS

相关帖子

相关日志

关闭 安全验证

标签: EMS

相关帖子

相关日志

关闭安全验证