用python做生信_1 python生信入门

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 格式 四行