1. 解决的实际问题,计算人类参考基因组hg38每条染色体碱基数据ATCG的比例
2. 计算人类参考基因组hg38的外显子区域碱基数量
3. 目的是了解基础python语法
4. 参考书目:python编程从入门到实践,参考视频资料见最底部
打开python3
vim 1.txt
:wq ##Esc+ :wq保存退出
python3 # 打开python3
f = open('./1.txt','r')
for i in f:
print (i)
with open('./1.txt') as file_object:
for line in file_object:
print (line)
图1
f = open('./2.txt','w')
with open('./1.txt','r') as imput_f:
for line in imput_f:
line = line.strip()
result= int(line)+int(line)
f.write(str(result)+'\n')
f.close()
图2
output_file = open('./2.txt','a') ## 附加内容
with open('./1.txt','r') as imput_f:
for line in imput_f:
line = line.strip()
result= int(line)+int(line)
output_file.write(str(result)+'\n')
output_file.close()
图3
完成探索人类基因组GC含量与gene密度的关系;
提示:
1.将染色体成分1Mb大小的bin统计GC含量;
2.将染色体分为1Mb大小的bin统计gene的密度;
3.使用R将这两者用散点图进行绘图;
4.使用R的线性回归函数lm()进行回归;
5.R回归部分,参考资料
blog.fens.me/r-linear-regression/
def my_abs(imput):
if imput >= 0:
return imput
else:
return -imput
my_abs(-10)
作业
人类的线粒体长度大约为16.7Kb,请读取人类线粒体参考基因组序列信息(fasta格式),统计各种碱基的数量以及人类线粒体参考基因组的总长度。
编程作业
# _*_ coding:UTF-8 _*_
### 保证中文不报错pycharm
选取参考基因组hg38作为测试数据1000-1300行 n 和 p是同时出现的
cat ~/reference/UCSC/hg38/hg38.fa| sed -n '1000,1300p' >> ~/project/python/t.txt
cd - # 回到上一个目录
python3
图4
chrM_genome = ''
with open('./t.txt','r') as imput_f:
for line in imput_f:
line = line.strip().upper() ## 去除前后空格,大写
chrM_genome = chrM_genome + line
print (len(chrM_genome))
print (chrM_genome.count('A'))
print ('the char A count is : {0}'.format(chrM_genome.count('A'))) # format 将其格式化
Fastq 格式 四行