1 python生信入門

1. 解決的實(shí)際問題，計(jì)算人類參考基因組hg38每條染色體堿基數(shù)據(jù)ATCG的比例

2. 計(jì)算人類參考基因組hg38的外顯子區(qū)域堿基數(shù)量

3. 目的是了解基礎(chǔ)python語法

4. 參考書目：python編程從入門到實(shí)踐，參考視頻資料見最底部

5. 測試數(shù)據(jù)下載地址鏈接：https://share.weiyun.com/5Xmmzdz 密碼：6q2y5j

打開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')  ## 附加內(nèi)容

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密度的關(guān)系；

提示：

1.將染色體成分1Mb大小的bin統(tǒng)計(jì)GC含量；

2.將染色體分為1Mb大小的bin統(tǒng)計(jì)gene的密度髓迎；

3.使用R將這兩者用散點(diǎn)圖進(jìn)行繪圖；

4.使用R的線性回歸函數(shù)lm()進(jìn)行回歸；

5.R回歸部分狮杨，參考資料

blog.fens.me/r-linear-regression/

https://github.com/menghaowei

def my_abs(imput):
     if imput >= 0:
         return imput
     else:
         return -imput
        
my_abs(-10)

作業(yè)

人類的線粒體長度大約為16.7Kb，請讀取人類線粒體參考基因組序列信息（fasta格式）到忽，統(tǒng)計(jì)各種堿基的數(shù)量以及人類線粒體參考基因組的總長度橄教。

編程作業(yè)

# _*_ coding:UTF-8 _*_
### 保證中文不報(bào)錯(cuò)pycharm


選取參考基因組hg38作為測試數(shù)據(jù)1000-1300行 n 和 p是同時(shí)出現(xiàn)的
cat ~/reference/UCSC/hg38/hg38.fa| sed -n '1000,1300p' >> ~/project/python/t.txt
cd - # 回到上一個(gè)目錄
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 格式 四行

圖5
>>> with open(r'./fa_q.txt','r') as imput:
...  for index,line in enumerate(imput):
...   if index % 4==0:
...    print ('>'+line.strip()[1:])
...   elif index % 4 ==1:
...    print (line.strip())
...   elif index % 4==2:
...    continue
...   elif index %4==3:
...    continue
... 
>SRR1039512.325 HWI-ST177:314:C12K6ACXX:1:1101:5657:2438 length=63
TCAGAAGAAAAGCATATGTCTGCCATGGGACTATTGCAGTGCGTCTCCATCAGTGTTAACACA
>SRR1039512.326 HWI-ST177:314:C12K6ACXX:1:1101:5550:2446 length=63
TAAGCAGTGCTTGAATTATTTGGTTTCGGTTGTTTTCTATTAGACTATGGTGAGCTCAGGTGA
>SRR1039512.327 HWI-ST177:314:C12K6ACXX:1:1101:5908:2367 length=63
TGGGGAATAGCTTCTCAAGGGCAAGTGCTGTTTCCATACTAGTTCTTTTCCTTGCTCTCTTAT
>SRR1039512.328 HWI-ST177:314:C12K6ACXX:1:1101:5852:2396 length=63
GGAGGAAGCCCACGGTGGGATCTGGCACAAAGATGTAGAGCCGTTTCCGCTCATCGGTCTCCA
>SRR1039512.329 HWI-ST177:314:C12K6ACXX:1:1101:5768:2409 length=63
ACCTGGCAGTGCTGCAACAGTATCTGCCTCTACAAGTAACATCATACCCCCAAGACACCAGAA
>SRR1039512.330 HWI-ST177:314:C12K6ACXX:1:1101:5840:2413 length=63
GCAGAGGTAATTGCCTGAATCCTTCTGTTGTAGACTACGTAGCAGAAGGCCTTGATCTGTCCT
>SRR1039512.331 HWI-ST177:314:C12K6ACXX:1:1101:5826:2463 length=63
GAAAGGCAAGGAGGAAGCTTATCTATGAAAAAGCAAAGCACTATCACAAGGAATATAGGC
>SRR1039512.332 HWI-ST177:314:C12K6ACXX:1:1101:5888:2477 length=63
GTCTGTCTCTTGCTTCAACAGTGTTTGGACGGAACAGATCCGGGGACTCTCTTCCAGCCTCCG
>SRR1039512.333 HWI-ST177:314:C12K6ACXX:1:1101:6208:2288 length=63
GTCATAACCACCACCGCTTACACCTGGAGGTCCAGGAGGGCCAGGGGGACCAGGGGGGCCAGC

##每隔40一行
with open(r'C:\Users\kexin\Desktop\python\fa_q.txt','r') as imput:
  for index,line in enumerate(imput):
   if index % 4==0:
    print ('>'+line.strip()[1:])
   elif index % 4 ==1:
    for i in range(0,len(line.strip()),40):
     print (line.strip()[i:i+40])
   elif index % 4==2:
    continue
   elif index %4==3:
    continue

output_file = open (r'C:\Users\kexin\Desktop\python\result2.fastq','w')
with open(r'C:\Users\kexin\Desktop\python\fa_q.txt','r') as imput:
  for index,line in enumerate(imput):
   if index % 4==0:
    output_file.write(('>'+line.strip()[1:])+'\n')
   elif index % 4 ==1:
    for i in range(0,len(line.strip()),40):
     output_file.write(line.strip()[i:i+40]+'\n')
   elif index % 4==2:
    print (index)
   elif index %4==3:
    continue
output_file.close()

讀取fasta文件

with open (r'C:\Users\kexin\Desktop\python\result2.fastq','r') as input_file:
    seq=''
    header= input_file.readline().strip()[1:]  ## 讀取第一行
    for line in input_file:
        if line[0]!='>':
            seq =seq.strip()+line  ### 這里是要確保剛剛被40個(gè)一行分隔的seq數(shù)據(jù)重新疊加
        elif line[0]=='>':
            print (header)
            print (seq)
            header = line.strip()[1:]
            seq = ''       ### 一個(gè)循環(huán)結(jié)束將seq清空
    print (header)
    print (seq)

一個(gè)能讀取fasta格式的函數(shù)（收藏版）


def read_fasta(file_path=''):
    line = ''
    try:
        fasta_handle=open(file_path,'r')
    except:
        raise IOError('you input FAST file is not right')
    while True:
        line = fasta_handle.readline()
        if line == '':
            return
        if line[0] =='>':
            break
 ## while the file is not empty, load FASTA file   
    while True:
        if line[0] !='>':
            raise ValueError("Records in Fasta files should start with '>'")
        title = line[1:].rstrip()
        lines =[]
        line = fasta_handle.readline()
        while True:
            if not line:
                break
            if line[0] == '>':
                break
            lines.append(line.rstrip())
            line = fasta_handle.readline()
        yield title, ''.join(lines).replace('','').replace('\r','')
        if not line:
            return
    fasta_handle.close()
    assert False, 'your input fasta file have format problem'
    
    
read_fasta(file_path=r'C:\Users\kexin\Desktop\python\result2.fastq')  ## 出現(xiàn)一個(gè)報(bào)錯(cuò)喘漏，在windows里面护蝶，要加上r，不然\代表轉(zhuǎn)義:參考https://blog.csdn.net/caibaoH/article/details/78335094


for header,seq in read_fasta(file_path=r'C:\Users\kexin\Desktop\python\result2.fastq'):
    print (header)
    print (seq)

計(jì)算人類參考基因組的序列長度

分析步驟：

利用讀取fasta的代碼讀取基因組序列信息

統(tǒng)計(jì)每條序列的長度翩迈，并輸出

統(tǒng)計(jì)人類參考基因組的總長度持灰，并輸出

統(tǒng)計(jì)每條染色體N的長度，并輸出

提取基因組特定區(qū)域的序列负饲，并輸出（支持+/-鏈）

hg38_genome={}
for chr_name, chr_seq in read_fasta(file_path=r'/four/yqchen/reference/UCSC/hg38/hg38.fa'):
    hg38_genome[chr_name]=chr_seq

這里要用絕對路徑不然報(bào)錯(cuò)

圖6
hg38_genome.keys()
for chr_name in hg38_genome.keys():
    print (chr_name,len(hg38_genome.get(chr_name)))
## 這里有點(diǎn)不對堤魁，但是無所謂了
圖7
hg38_total_len=0
for index in range (1,26):
    if index<= 22:
        chr_name = 'chr{0}'.format(index)
    elif index == 23:
        chr_name = 'chrX'
    elif index == 24:
        chr_name = 'chrY'
    elif index == 25:
        chr_name = 'chrM'
    print (chr_name,len(hg38_genome.get(chr_name)))

這里結(jié)果是正確的，這是為什么返十？這里沒搞明白的妥泉。

圖8
hg38_total_len = 0
for index in range (1,26):
    if index<= 22:
        chr_name = 'chr{0}'.format(index)
    elif index == 23:
        chr_name = 'chrX'
    elif index == 24:
        chr_name = 'chrY'
    elif index == 25:
        chr_name = 'chrM'
    hg38_total_len= hg38_total_len + int(len(hg38_genome.get(chr_name))) 

print (hg38_total_len) ### 和上面那一句分開執(zhí)行
圖9

圖10
len(hg38_genome['chr21'])
hg38_genome['chr21'].count('N')
hg38_genome['chr22'].count('N')

## 正負(fù)鏈

chr_name='chr1'
chr_start=10000
chr_end=10300
hg38_genome[chr_name][chr_start-1:chr_end-1].upper()


#### 反向互補(bǔ)
import string
translate_table=string.maketrans('ATGC','TACG')
a=hg38_genome[chr_name][chr_start-1:chr_end-1].upper()
a.translate(translate_table)[::-1]


##方向互補(bǔ)定義函數(shù)
import string
def com_rev(input):
    translate_table=string.maketrans('ATGC','TACG')
    return input.upper().translate(translate_table)[::-1]

編碼基因長度，從參考基因組注釋信息找

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()
    for line in input_f:
        line=line.strip().split('\t')
        print (line)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        exon_count = int(line [8])
        exon_start =line[9].strip(',') ### 去掉最后的逗號
        print (exon_start)
        break



with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start =line[9].strip(',')
        exon_end =line[10].strip(',')### 去掉最后的逗號
        print (exon_start)
        print (exon_end)
        break


with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start = line[9].strip(',').split(',') ### map函數(shù)int用于后面的列表
        exon_end = line[10].strip(',').split(',')
        print (exon_start)
        print (exon_end)
        break




with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_count = int(line [8])
        exon_start = map(int,line[9].strip(',').split(',')) ### map函數(shù)int用于后面的列表
        exon_end = map(int,line[10].strip(',').split(','))
        print (exon_start)
        print (exon_end)
        break

with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = map(int,line[9].strip(',').split(',')) ### map函數(shù)int用于后面的列表
        exon_end = map(int,line[10].strip(',').split(','))
        
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        print (exon_total_len)
        print (exon_start)
        print (exon_end)
        break

這里python3 會(huì)報(bào)錯(cuò)

1555423217495.png

解決方案洞坑，在map基礎(chǔ)上加上list

https://blog.csdn.net/mingyuli/article/details/81238858

1555423298540.png
with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = list(map(int, line[9].strip(',').split(','))) ### map函數(shù)int用于后面的列表
        exon_end = list(map(int, line[10].strip(',').split(',')))
        
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        print (exon_total_len)
        print (exon_start)
        print (exon_end)
        break

完整計(jì)算版：

gene_exon_len_dict={}
with open(r'C:/Users/kexin/Desktop/python/hg38_RefSeq_genetable') as input_f:
    header = input_f.readline()   #### 讀一行 不要用readlines
    for line in input_f:
        line=line.strip().split('\t')
        # print (line)
        exon_total_len=0
        exon_count = int(line [8])
        exon_start = list(map(int, line[9].strip(',').split(','))) ### map函數(shù)int用于后面的列表
        exon_end = list(map(int, line[10].strip(',').split(',')))
         
        for index in range(exon_count):
            exon_total_len= exon_total_len+exon_end[index]-exon_start[index]
        gene_id = line[12]
        
        if gene_exon_len_dict.get(gene_id) is None: #如果gene_exon_len_dict字典里沒有這個(gè)基因盲链，加入這個(gè)基因的長度
            gene_exon_len_dict[gene_id] = exon_total_len
        else:
            if gene_exon_len_dict.get(gene_id) < exon_total_len:
                gene_exon_len_dict[gene_id] = exon_total_len   


len(gene_exon_len_dict)

exon_total =0 
for gene_id in gene_exon_len_dict.keys():
    exon_total = exon_total+gene_exon_len_dict[gene_id]
    # exon_total = exon_total+gene_exon_len_dict.values(gene_id) 這是錯(cuò)誤寫法
    
print (exon_total)


exon_rate = 111336157 / 3088286401 *100
exon_rate

推薦視頻

https://www.bilibili.com/video/av16064125?from=search&seid=5895915127692784448

https://www.bilibili.com/video/av15887563?from=search&seid=5895915127692784448

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子刽沾，更是在濱河造成了極大的恐慌本慕，老刑警劉巖，帶你破解...
    沈念sama閱讀 218,525評論 6贊 507
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件悠轩，死亡現(xiàn)場離奇詭異间狂，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)火架，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,203評論 3贊 395
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門鉴象，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人何鸡，你說我怎么就攤上這事纺弊。” “怎么了骡男？”我有些...
    開封第一講書人閱讀 164,862評論 0贊 354
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵淆游，是天一觀的道長。 經(jīng)常有香客問我隔盛，道長犹菱，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 58,728評論 1贊 294
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任吮炕，我火速辦了婚禮腊脱，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘龙亲。我一直安慰自己陕凹，他們只是感情好，可當(dāng)我...
    茶點(diǎn)故事閱讀 67,743評論 6贊 392
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布鳄炉。 她就那樣靜靜地躺著杜耙，像睡著了一般。 火紅的嫁衣襯著肌膚如雪拂盯。 梳的紋絲不亂的頭發(fā)上佑女，一...
    開封第一講書人閱讀 51,590評論 1贊 305
  • 城市分裂傳說
    那天，我揣著相機(jī)與錄音磕仅，去河邊找鬼珊豹。 笑死，一個(gè)胖子當(dāng)著我的面吹牛榕订，可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播蜕便，決...
    沈念sama閱讀 40,330評論 3贊 418
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼劫恒，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起两嘴，我...
    開封第一講書人閱讀 39,244評論 0贊 276
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤丛楚，失蹤者是張志新（化名）和其女友劉穎，沒想到半個(gè)月后憔辫，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體趣些，經(jīng)...
    沈念sama閱讀 45,693評論 1贊 314
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,885評論 3贊 336
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年贰您，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了坏平。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 40,001評論 1贊 348
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡锦亦，死狀恐怖舶替，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情杠园，我是刑警寧澤顾瞪，帶...
    沈念sama閱讀 35,723評論 5贊 346
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站抛蚁，受9級特大地震影響陈醒，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜瞧甩，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,343評論 3贊 330
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一钉跷、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧亲配，春花似錦尘应、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,919評論 0贊 22
  • 一樁弒父案犬钢，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至思灰，卻和暖如春玷犹，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背洒疚。 一陣腳步聲響...
    開封第一講書人閱讀 33,042評論 1贊 270
  • 情欲美人皮
    我被黑心中介騙來泰國打工歹颓， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人油湖。 一個(gè)月前我還...
    沈念sama閱讀 48,191評論 3贊 370
  • 代替公主和親
    正文 我出身青樓巍扛，卻偏偏與公主長得像，于是被迫代替她去往敵國和親乏德。 傳聞我的和親對象是個(gè)殘疾皇子撤奸，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,955評論 2贊 355