前言
比對軟件很多艳狐,首先大家去收集一下聘殖,因為我們是帶大家入門,請統(tǒng)一用hisat2续挟,并且搞懂它的用法紧卒。直接去hisat2的主頁下載index文件即可,然后把fastq格式的reads比對上去得到sam文件诗祸。接著用samtools把它轉(zhuǎn)為bam文件常侦,并且排序(注意N和P兩種排序區(qū)別)索引好,載入IGV贬媒,再截圖幾個基因看看聋亡!順便對bam文件進(jìn)行簡單QC,參考直播我的基因組系列际乘。生信技能樹
參考基因組及注釋
由于在第一節(jié)已經(jīng)準(zhǔn)備好了相應(yīng)的軟件坡倔,只需要再準(zhǔn)備好相應(yīng)的數(shù)據(jù)就行了。在HISAT2官網(wǎng)即可下載現(xiàn)成的index文件脖含,這里使用小鼠常用的mm10罪塔。
接下來就用下載相應(yīng)的注釋,切記不同及基因組版本的對應(yīng)關(guān)系养葵。
通過gencode下載最新小鼠GRCm38版本的gtf注釋征堪。這個頁面也收集好了不同分類的基因集合,比如lncRNA的注釋关拒。
序列比對
ref=/mnt/e/0ngs/ref/mm10_hisat_index/mm10 #index
# 根據(jù)注釋文件提取已知splices位點信息
hisat2_extract_splice_sites.py gencode.vM13.annotation.gtf > splicesites.txt
# 開始比對
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
比對結(jié)果統(tǒng)計
其中佃蚜,concordantly指方向和相對距離都與建庫一致;disconcordantly指配對的兩個read都能比對上着绊,但方向或者(和)相對距離不符合谐算。s4、s5指不能同時比對上的配對read归露,拆分按單端測序條件比對的情況洲脂。
Sort
### sort -n(read name排序) ###
# name排序方便提高后續(xù)read count的效率
samtools sort -n -@ 6 Akap95_1.sam -o Akap95_1.Nsort.bam
# position排序利于存儲壓縮
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