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

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

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

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

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

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

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

言归正传，今天我们要继续来聊一聊上次未完的话题，小麦穗粒数转录组分析(一),上次我们说到筛选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千个位点可以拿来做分析了。今天就先到这，下周我们再约。