数据下载来自ncbi的ftp，想做这个的原因来自于jimmy的生信技能树的论坛帖子：
http://www.biotrainee.com/thread-265-1-1.html
我们都知道，一个人的全基因组测序会与参考基因组hg19/hg38有300~600万的snp差异，那么现在已有测序的所有样本与参考基因组hg19/hg38不同的地方是多少呢？全基因组测序已经超过10万例了，对snp位点记录最详细的当属 dbsnp啦，很容易去NCBI的ftp里面下载最新版的dbsnp文件，约15G大小，最新版的话！
如果只需要看统计信息，其实在dbsnp的官网主页即可看到。但是我们这个数据处理任务，所以必须去下载该文件，对文件进行解析。除了可以看有多少位点，还可以探索更多！欢迎爱学习的小伙伴在下面跟帖讨论自己的学习结果。
1.dbsnp数据库的下载
我下载的是下面这个地址的vcf.gz格式的，实在找不到txt格式的在哪，于是就下这个了。压缩文件大小是7个多G，解压之后有50多G，估计是包含了很多的注释信息。也没有选择对应于哪个参考基因组的snp，估计这个是最全的。对于这样的大文件，我习惯于IBM的ascp下载了(当然，wgetg一样可以)，代码如下：
ascp -QT -l 6M -k1 -i asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.nih.gov:/snp/organisms/human9606/VCF/All20170710.vcf.gz D:/
2.数据库文件的处理和snp的计数
shell脚本应该也是可以搞定的，我对于shell是门外汉，先用python解决问题，尽管代码要多好几倍，估计shell也就一行的事。哪天学会了再补上。
处理思路是，发现文件里除了#开头注释行，其余均为数据，所以排除#开头的行，然后把染色体号，rs号，正常基因和变异基因提取到一个新的文件，代码如下：
fin = open('D:\\All_20170710.vcf')
fout = open('D:\\all_snp.csv', 'w')
i =0
for line in fin:
if line[0] !='#':
t = line.strip().split('\t')
line = t[0] +'\t'+ t[2] + '\t' + t[3] + '\t' + t[4]
i += 1
#print(i, t)
#print(str(i) + '\t' + line + '\n')
fout.write(str(i) + '\t' + line + '\n')
#break
print(i)
fout.close()
#运行结果如下：
325174796
[Finished in 2792.9s]
可以看出，总数竟然高达325，174，796，这是三亿多呀，远远超出了我的想象，我以为只会有千万以内呢！但是，问题是现在能确定功能的snp位点只有上千个左右，也就是基因检测报告应用的数目，科技仍需加油呀！
3.统计各染色体的分布情况
用个字典解决问题吧！下面是是我的代码：
fin = open('D:\\all_snp.csv')
di = {}
chromesome = '1'
no = 0
for line in fin:
#print(line.strip().split('\t')[1])
#print(chromesome)
if line.strip().split('\t')[1] == chromesome:
no += 1
#print(no)
else:
di[chromesome] = no
no = 0
chromesome = line.strip()[1]
#break
print(di)
#运行结果
{'1': 25140135, '2': 27537299, '3': 22606013, '4': 21849502, '5': 20438556, '6': 19276947, '7': 17929935, '8': 17564397, '9': 13893120, '10': 15064971, '11': 15437930, '12': 15007421, '13': 10958764, '14': 10218626, '15': 9256140, '16': 10351834, '17': 9133768, '18': 8643026, '19': 7179582, '20': 7100293, '21': 4229709, '22': 4252450, 'X': 11733164, 'Y': 368898}
[Finished in 537.0s]
这个运行时间还算可以接受哈，500多秒，也就是十分钟不到。那么，排个序，哪个染色体最多，我猜应该和染色体大小一致。
3.多少排序
chromesome = {'1': 25140135, '2': 27537299, '3': 22606013, '4': 21849502, '5': 20438556, '6': 19276947, '7': 17929935, '8': 17564397, '9': 13893120, '10': 15064971, '11': 15437930, '12': 15007421, '13': 10958764, '14': 10218626, '15': 9256140, '16': 10351834, '17': 9133768, '18': 8643026, '19': 7179582, '20': 7100293, '21': 4229709, '22': 4252450, 'X': 11733164, 'Y': 368898}
li = sorted(chromesome.items(), key= lambda item:item[1], reverse = True)
print(li)
#运行结果
[('2', 27537299), ('1', 25140135), ('3', 22606013), ('4', 21849502), ('5', 20438556), ('6', 19276947), ('7', 17929935), ('8', 17564397), ('11', 15437930), ('10', 15064971), ('12', 15007421), ('9', 13893120), ('X', 11733164), ('13', 10958764), ('16', 10351834), ('14', 10218626), ('15', 9256140), ('17', 9133768), ('18', 8643026), ('19', 7179582), ('20', 7100293), ('22', 4252450), ('21', 4229709), ('Y', 368898)]
[Finished in 0.2s]
可以看出2号染色体最多，那么看一下染色体的大小情况，来自论坛编程直播课程作业贴：
chr GC_ratio N_ratio Length N A T C G
>chr1 0.41744 0.09617 249250621 23970000 65570891 65668756 47024412 47016562
>chr2 0.40244 0.02054 243199373 4994855 71102632 71239379 47915465 47947042
>chr3 0.39694 0.01629 198022430 3225295 58713343 58760485 38653197 38670110
>chr4 0.38248 0.01827 191154276 3492600 57932980 57952068 35885806 35890822
>chr5 0.39516 0.0178 180915260 3220000 53672554 53804137 35089383 35129186
>chr6 0.39611 0.02174 171115067 3720001 50554433 50533923 33143287 33163423
>chr7 0.40751 0.02378 159138663 3785000 45997757 46047257 31671670 31636979
>chr8 0.40176 0.02374 146364022 3475100 42767293 42715025 28703983 28702621
>chr9 0.41317 0.14921 141213431 21070000 35260078 35243882 24826212 24813259
>chr10 0.41585 0.03114 135534747 4220009 38330752 38376915 27308648 27298423
>chr11 0.41566 0.02872 135006516 3877000 38307244 38317436 27236798 27268038
>chr12 0.40812 0.02518 133851895 3370502 38604831 38624517 26634995 26617050
>chr13 0.38527 0.17001 115169878 19580000 29336945 29425459 18412698 18414776
>chr14 0.40887 0.17755 107349540 19060000 25992966 26197495 18027132 18071947
>chr15 0.42201 0.20322 102531392 20836626 23620876 23597921 17247582 17228387
>chr16 0.44789 0.12694 90354753 11470000 21724083 21828642 17630040 17701988
>chr17 0.4554 0.04187 81195210 3400000 21159933 21206981 17727956 17700340
>chr18 0.39785 0.0438 78077248 3420019 22465380 22489493 14838685 14863671
>chr19 0.4836 0.05615 59128983 3320000 14390632 14428951 13478255 13511145
>chr20 0.44126 0.05585 63025520 3520000 16523053 16725227 13107828 13149412
>chr21 0.40833 0.27059 48129895 13023253 10422924 10348785 7160212 7174721
>chr22 0.47988 0.31985 51304566 16410021 9094775 9054551 8375984 8369235
>chrX 0.39496 0.02686 155270560 4170000 45648952 45772424 29813353 29865831
>chrY 0.39965 0.56793 59373566 33720000 7667625 7733482 5099171 5153288
发现规律并不完全一致。好了，探索先到这里。
