2017-07-04

今天开始跟着公司的一个人做流程了，他也是第一次做。流程是先做一遍正常的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模式，百度去百度去、、、、、、、、、、、不过，成功添加代码块，啦啦啦啦啦。。