今天開(kāi)始跟著公司的一個(gè)人做流程了胧华，他也是第一次做温数。流程是先做一遍正常的RNA-seq的數(shù)據(jù)。
序列比對(duì)使用的軟件是HISAT2軟件驾荣，和bowtie2是一樣的外构。
第一步是下載數(shù)據(jù)，ensemble上下載的人類(lèi)基因組fasta文件
第二步播掷，使用hisat2-bulid把基因組切分成小的index审编，這樣比對(duì)的時(shí)候能提高效率
命令行:hisat2-bulid genome.fa index
之后會(huì)得到一個(gè)index的文件夾
第三步，使用hisat把雙端測(cè)序數(shù)據(jù)比對(duì)到基因組上
命令行:hisat2 -p 8 --dta -x index -1 第一端測(cè)序數(shù)據(jù).fastq -2 第二端測(cè)序數(shù)據(jù).fastq -S 輸出結(jié)果文件.sam
參數(shù)詳解http://blog.sciencenet.cn/blog-759995-990471.html
大概看一下sam文件歧匈，一般有用的分別是第二個(gè)值垒酬，代表正負(fù)鏈或者沒(méi)匹配之類(lèi)的；第三個(gè)值染色體件炉，正常是1到22外加X(jué)Y勘究，像如圖這種奇形怪狀的就是沒(méi)匹配的；第四個(gè)值斟冕，代表位置信息口糕。
第四步，畫(huà)圖磕蛇，對(duì)每天染色體景描，每100k一個(gè)區(qū)間，計(jì)算map到每個(gè)區(qū)間上的read數(shù)秀撇，正負(fù)鏈畫(huà)一張圖上
因?yàn)閿?shù)據(jù)沒(méi)給我超棺，我就自己瞎隨機(jī)了正負(fù)鏈，哈哈(ω)
import pandas as pd
import matplotlib.pyplot as plt
import random as rd
c = [str(x) for x in range(1,23)] #這里不用str下面賦值1-22都是空
c.append('X')
c.append('Y')
data = pd.DataFrame(None,index = c,columns = ['forward','reverse'])
for x in data.index:
data.loc[x,'forward'] = []
data.loc[x,'reverse'] = []
f = open("E:/geneX/0704/index-TAGCTT_H35YMCCXY_L3_tout_accepted_hit.sam")
for line in f.readlines():
if(line[0] == '@'):
next
else:
l = line.split('\t') #l[2]:chr l[3]:position
if(l[2] in c):
if rd.random()>0.5:
data.loc[l[2],'forward'].append(int(int(l[3])/100000.))
else:
data.loc[l[2],'reverse'].append(int(int(l[3])/100000.))
這樣得到的data中index是染色體1-22呵燕、X棠绘、Y，columns是forward和reverse再扭，只不過(guò)是我隨機(jī)噠氧苍。
之后，可以畫(huà)圖啦~~~~
for x in data.index:
forward = pd.DataFrame(data.loc[x]['forward'],columns = ['posi']).groupby('posi').size()
reverse = pd.DataFrame(data.loc[x]['reverse'],columns = ['posi']).groupby('posi').size()
x = forward.index
y = forward.values
plt.plot(x,y)
x2 = reverse.index
y2 = reverse.values
plt.plot(x2,-y2)
plt.show()
break;
看一張圖泛范，不得不說(shuō)python的顏色還是很好看的候引，默認(rèn)參數(shù)就挺好看了
因?yàn)槭请S機(jī)的，所以正負(fù)鏈比較對(duì)稱(chēng)敦跌，正常不會(huì)出現(xiàn)對(duì)稱(chēng)情況的。
好啦明天可以去看個(gè)lncRNA的測(cè)序數(shù)據(jù)啦
我添真的是很棒棒，柠傍，麸俘，，狠狠棒棒惧笛，比如這個(gè)markdown模式从媚，百度去百度去、患整、拜效、、各谚、紧憾、、昌渤、赴穗、、膀息、不過(guò)般眉，成功添加代碼塊，啦啦啦啦啦潜支。甸赃。
