今天开始跟着公司的一个人做流程了，他也是第一次做。流程是先做一遍正常的RNA-seq的数据。
序列比对使用的软件是HISAT2软件，和bowtie2是一样的。
第一步是下载数据，ensemble上下载的人类基因组fasta文件
mmexport1499176610997.jpg
第二步，使用hisat2-bulid把基因组切分成小的index，这样比对的时候能提高效率
命令行:hisat2-bulid genome.fa index
之后会得到一个index的文件夹
第三步，使用hisat把双端测序数据比对到基因组上
命令行:hisat2 -p 8 --dta -x index -1 第一端测序数据.fastq -2 第二端测序数据.fastq -S 输出结果文件.sam
参数详解http://blog.sciencenet.cn/blog-759995-990471.html
大概看一下sam文件，一般有用的分别是第二个值，代表正负链或者没匹配之类的；第三个值染色体，正常是1到22外加XY，像如图这种奇形怪状的就是没匹配的；第四个值，代表位置信息。
第四步，画图，对每天染色体，每100k一个区间，计算map到每个区间上的read数，正负链画一张图上
因为数据没给我，我就自己瞎随机了正负链，哈哈(ω)
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，只不过是我随机哒。
之后，可以画图啦~~~~
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;
看一张图，不得不说python的颜色还是很好看的，默认参数就挺好看了
Paste_Image.png
因为是随机的，所以正负链比较对称，正常不会出现对称情况的。
好啦明天可以去看个lncRNA的测序数据啦
我添真的是很棒棒，，，，狠狠棒棒，比如这个markdown模式，百度去百度去、、、、、、、、、、、不过，成功添加代码块，啦啦啦啦啦。。
