生信指令碼練習（6） 求read每個位點的cg分布

這是另外一題，也是cg含量，

但是是求乙個fastq檔案中每個read位點的cg分布情況

這個解法是逐行讀取檔案，比較優雅

number  = {}

buffer = 200
for i in range(buffer):
kkk = i              
number[kkk] = 0
with
open("test1.fastq","r") as f:
line_number = 0
while true:       
line = f.readline()
line_number += 1
if line_number % 4 == 2 :             
for i in range(len(line)):
kkk = i 
ifline[kkk] == "c"
orline[kkk] == "g" :
number[kkk] = 1 + number[kkk] 
else:
pass
ifnotlen(line):
break
#print(number) 
f = open('cg.txt', "w")
f.write('location'+ "\t" + 'number'+ '\n' )
for k,v in
number.items():    
f.write(str(k+1)+"\t"+str(v)+"\n")    
f.close()

setwd("\\desktop") 

data
= t) 
plot(data
$location,data
$number,main=
"折線圖",xlab=
"number",ylab=
"len",col=
"red",cex=
1,lwd=
1,type
='o')

# 每條reads的cg分布 畫圖

library(ggplot2)
setwd("c:\\desktop") 
dataa 
#print(data)
jpeg("pie.jpeg",width=2000,height=2000,res=300) 
ggplot(dataa, aes(x = location, y = number)) +
geom_line()+
geom_point()
dev.off()

# 8.14 修改 這也是逐行讀，修正了突然出現多幾個鹼基的read就會報錯的bug

#另外用matplotlib畫了個圖,最後的轉折是因為沒去掉末尾的換行符。。
import os
import matplotlib.pyplot as plt
os.chdir("d:/") 
with
open("a1_r1.fq","r") as f:
number  = {}
line_number = 0
while true:       
line = f.readline().strip()
line_number += 1
if line_number % 4 == 2 :             
for i in range(len(line)):
if i in
number:
ifline[i] == "c"
orline[i] == "g" :
number[i] = 1 + number[i]    
else:
ifline[i] == "c"
orline[i] == "g" :
number[i] = 1
else:
number[i] = 0
iflen(line) == 0:
break
print(number)
line_number = (line_number-1)/4 
fanal = {}
for k,v in
number.items():  
fanal[k+1] = v/line_number 
x = list(fanal.keys())
y = list(fanal.values())
plt.plot(x, y)
plt.ylabel('frequency')
plt.xlabel('location')
plt.title('cg of reads')
#plt.text(6, .15, r'$\mu=100,\ \sigma=15$')
plt.show()

生信基礎（二） 生信學習資料

原創 hxj7 上次談到生信人員需要熟練掌握一些程式語言，還講了perl和python的選擇問題。那麼，如果已經選定了一門程式語言，到底該如何學習它呢？今天的我們可以通過mooc跟著名師學習或者上知乎提問，幸運的話還能得到大牛指點。不過，在我剛接觸程式設計的時候，mooc和知乎都還未興起，所以我都是...

常見生信操作

安裝samtools conda install samtools srand 隨機數發生器。設定固定的種子，保證每次出來的結果一致 rand 返回 0,1 之間的隨機數，包含0不包含1 1.產生隨機的基因組檔案 echo 1 awk v seed 1 v label chr v chrnum 4 ...

生信轉崗心得

應健明的邀請，我也寫一篇關於轉行做生物資訊的心得，本來以為很輕鬆就可以寫出來的，但是發現並不那麼好寫。如果放在去年剛轉崗之時，我想應該更順手，那時感觸良多且深刻。不過我雖是個健忘的人，但是還清醒地記得去年7月份通過公司的轉崗答辯之時，心情無比的愉快與美麗，覺得終於從一名生信愛好者成為一名正規軍。之前...