前言

比對(duì)軟件很多，首先大家去收集一下，因?yàn)槲覀兪菐Т蠹胰腴T(mén)，請(qǐng)統(tǒng)一用hisat2，并且搞懂它的用法。直接去hisat2的主頁(yè)下載index文件即可，然后把fastq格式的reads比對(duì)上去得到sam文件。接著用samtools把它轉(zhuǎn)為bam文件，并且排序(注意N和P兩種排序區(qū)別)索引好，載入IGV，再截圖幾個(gè)基因看看！順便對(duì)bam文件進(jìn)行簡(jiǎn)單QC，參考直播我的基因組系列。生信技能樹(shù)

參考基因組及注釋

由于在第一節(jié)已經(jīng)準(zhǔn)備好了相應(yīng)的軟件，只需要再準(zhǔn)備好相應(yīng)的數(shù)據(jù)就行了。在HISAT2官網(wǎng)即可下載現(xiàn)成的index文件，這里使用小鼠常用的mm10。

1503648924351.png

接下來(lái)就用下載相應(yīng)的注釋，切記不同及基因組版本的對(duì)應(yīng)關(guān)系。

1503649711228.png

1503649797663.png

通過(guò)gencode下載最新小鼠GRCm38版本的gtf注釋。這個(gè)頁(yè)面也收集好了不同分類的基因集合，比如lncRNA的注釋。

序列比對(duì)

ref=/mnt/e/0ngs/ref/mm10_hisat_index/mm10 #index
# 根據(jù)注釋文件提取已知splices位點(diǎn)信息
hisat2_extract_splice_sites.py gencode.vM13.annotation.gtf > splicesites.txt
# 開(kāi)始比對(duì)
for i in `seq 59 62`; do hisat2 -p 6 -x $ref --known-splicesite-infile ref/splicesites.txt -1 SRR35899$i_1.fastq.gz -2 SRR35899$i_2.fastq.gz -S SRR35899$i.sam 2>SRR35899$i.log; done 

比對(duì)結(jié)果統(tǒng)計(jì)

1503912458951.png

image.png

其中，concordantly指方向和相對(duì)距離都與建庫(kù)一致；disconcordantly指配對(duì)的兩個(gè)read都能比對(duì)上，但方向或者(和)相對(duì)距離不符合。s4、s5指不能同時(shí)比對(duì)上的配對(duì)read，拆分按單端測(cè)序條件比對(duì)的情況。

Sort

###  sort -n(read name排序)  ###
# name排序方便提高后續(xù)read count的效率
samtools sort -n -@ 6 Akap95_1.sam -o Akap95_1.Nsort.bam
# position排序利于存儲(chǔ)壓縮
samtools sort -@ 6 Akap95_1.sam -o Akap95_1.Psort.bam   #染色體坐標(biāo)排序
# index
ls *.Psort.bam |while read id;do (samtools index -@ 4 $id &);done