mashengwei的个人博客分享 http://blog.sciencenet.cn/u/mashengwei

博文

小麦穗粒数转录组分析(二)——SNP筛选

已有 6135 次阅读 2018-2-27 14:49 |系统分类:科研笔记| 小麦, SNP筛选

2/28 本期作者:Neal

小麦穗粒数转录组分析(二)

胖丫总爱胡思乱想,今天又问我“科研的意义是什么?”。

这个问题我可以回答,不过你要回答我另外一个问题,活着的意义是什么?

听到这个,胖丫已经无语了,故意用中指推了推鼻梁上的眼镜,或许这就是无声的抗议吧。

其实,我还想再问一个问题,意义的意义又是什么?

地球上生命多彩斑斓,也许只有我们人类会追问这类问题。也许活着的意义就是活着本身。

胖丫若有所思又茫茫然的点了下头,一言不发的转身去忙了。

言归正传,今天我们要继续来聊一聊上次未完的话题,小麦穗粒数转录组分析(一),上次我们说到筛选SNP了。

下面这一步就是筛选SNP了,这里使用的是SnpSift

cat 90_mini_core_UG.vcf | java -jar /data1/snpEff_latest_core/snpEff/SnpSift.jar filter "(QUAL > 30) &(MQ > 40) & (QD > 5) & (FS < 30.0) & (AN > 90) & (DP> 90) &(countVariant() > 9) & (countRef() > 9)"                 > 90_mini_core_UG_first_filter.vcf

使用下面的脚本再次筛选,筛选标准是,首先将杂合数据处理成missing data,单个品种内,DP小于3也处理成missing data,最后只保留具有2个等位变异的位点。这一步筛选之后还有大概33万个位点.

#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ ='wheatOmics'__author_email__ ='13148474750@163.com'with open('90_mini_core_UG_first_filter.vcf','r')as f:    for line in f:        if line.startswith('#'):            print line,        else:            lin = line.strip().split('')            if len(lin[4].split(','))==1:                for i in range(9,99):                    gt = lin[i].split(':')[0]                    if gt !='./.':                        ad1 = int(lin[i].split(':')[1].split(',')[0])                        ad2 = int(lin[i].split(':')[1].split(',')[1])                        dp = lin[i].split(':')[2]                        gq = lin[i].split(':')[3]                        pl = lin[i].split(':')[4]                        if int(dp)<5or max([ad1, ad2])/float(dp)<0.9:                            lin[i]='./.'+ lin[i][3:]                    if gt =='0/1':                        lin[i]='./.'+ lin[i][3:]                for i in lin:                    print i +'',                print'',

在上面的基础之上需要删除missing data大于20%的位点,同时MAF小于10%的位点也忽略。这一步筛选之后,就有点悲催了,只有5万5个位点用来进行下面的分析了。

#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ ='wheatOmics'__author_email__ ='13148474750@163.com'with open('90_mini_core_UG_first_second.vcf','r')as f:    for line in f:        if line.startswith('#'):            print line,        else:            lin = line.strip().split('')            missing =0            ref =0            alt =0            for i in range(9,99):                gt = lin[i].split(':')[0]                if gt =='./.':                    missing +=1                if gt =='0/0':                    ref +=1                if gt =='1/1':                    alt +=1            total = missing + ref + alt            min_value = min([ref, alt])            if total ==90and missing/float(total)<0.2and min_value/float(total)>=0.1:                for i in lin:                    print i +'',                print'',

下面就是染色体合并在一起。由于小麦每条染色体巨大,有些分析软件不支持,只能暂时分析时拆开,然后再合起来。

#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ ='wheatOmics'__author_email__ ='13148474750@163.com'chr =[['chr1A',471304005],['chr1B',438720154],['chr1D',452179604],['chr2A',462376173],['chr2B',453218924],       ['chr2D',462216879],['chr3A',454103970],['chr3B',448155269],['chr3D',476235359],['chr4A',452555092],       ['chr4B',451014251],['chr4D',451004620],['chr5A',453230519],['chr5B',451372872],['chr5D',451901030],       ['chr6A',452440856],['chr6B',452077197],['chr6D',450509124],['chr7A',450046986],['chr7B',453822637],       ['chr7D',453812268]]with open('IWGSCv1.0 90_mini_core_UG_first_third_filter.vcf','r')as f:    for line in f:        if line.startswith('#'):            print line,        else:            line = line.replace('_part1','')            line = line.strip().split('')            if line[0].endswith('part2'):                for i in chr:                    if line[0].split('_')[0]== i[0]:                        line[1]= int(line[1])+ int(i[1])                        line[0]= line[0].split('_')[0]            for m in line[:-1]:                print str(m)+'',            print line[-1]+'',

至于,SNP的功能功能注释,前期我们已经介绍过,使用SnpEff 对SNP结果进行分析。至此,这5万5千个位点可以拿来做分析了。今天就先到这,下周我们再约。

111322tqabi88lp9z87ap1.jpg




https://wap.sciencenet.cn/blog-1094241-1101453.html

上一篇:利用BSR-Seq技术快速定位小麦抗条锈病基因YrMM58和YrHY1
下一篇:小麦花药发育的形态学观察及发育时期判定
收藏 IP: 58.213.93.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-12-26 21:58

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部