一蓝谨、hisat2進(jìn)行轉(zhuǎn)錄組比對
1. 首先下載參考基因組
在NCBI,UCSC和Ensemble數(shù)據(jù)庫都可以下載到參考基因組膏燃,這里下載UCSC的參考基因組茂卦。
https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/
1.png
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
tar -zxvf chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*.fa #刪除合并前的文件
2. 參考基因組注釋下載
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz
gunzip gencode.v40.annotation.gtf.gz
2.png
3. 建立索引
將外顯子加入索引
hisat2_extract_exons.py gencode.v40.annotation.gtf > exons.txt
將剪切位點加入索引
hisat2_extract_splice_sites.py gencode.v40.annotation.gtf > ss.txt
建立索引
hisat2-build -p 4 --ss ss.txt --exon exons.txt hg19.fa hg19
注意:建立索引需要有較高的電腦配置或需要服務(wù)器,普通電腦一般去hisat2官網(wǎng)下載即可
下載地址: http://daehwankimlab.github.io/hisat2/download/#h-sapiens
3.png
4. 比對
hisat2 -p 4 -x /mnt/e/reference/hg19/genome -1 /mnt/d/bioinfo/project/rna/02clean/SRR15859344_1.clean.fastq.gz -2 /mnt/d/bioinfo/project/rna/02clean/SRR15859344_2.clean.fastq.gz -S /mnt/d/bioinfo/project/rna/03alin/SRR15859344.sam
p: 線程數(shù)
x: 索引位置
-1 -2: 雙端測序的兩端數(shù)據(jù)
-S: 輸出文件路徑
將sam文件轉(zhuǎn)為bam文件(samtools軟件)
單個文件處理
samtools sort -O bam -@ 5 -o SRR15859344.bam SRR15859344.sam
多個文件批量處理
ls *.sam |while read id; do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done
#".sam"表示去掉文件名后面的.sam组哩;
#basename $id ".sam".bam為去掉.sam加上.bam等龙,也就是把文件名后綴改為.bam处渣;
為bam文件構(gòu)建索引
為bam文件構(gòu)建索引
目的
必須對bam文件進(jìn)行默認(rèn)情況下的排序后,才能進(jìn)行index而咆。否則會報錯霍比。建立索引后將產(chǎn)生后綴為.bai的文件,用于快速的隨機處理暴备。很多情況下需要有bai文件的存在悠瞬,特別是顯示序列比對情況下。比如samtool的tview命令就需要涯捻;gbrowse2顯示reads的比對圖形的時候也需要浅妆。IGV顯示比對情況也需要。
單個文件處理
samtools index -@ 4 SRR15859344.bam
多個文件批量處理
ls *.bam |xargs -i samtools index {}
二障癌、counts表達(dá)矩陣
對bam文件進(jìn)行定量(featureCounts)
featureCounts -a /mnt/e/reference/gencode.v40.annotation.gtf -F gtf -t exon -g gene_id -p -T 4 -o /mnt/d/bioinfo/project/rna/04count/expr_matrix.txt /mnt/d/bioinfo/project/rna/03alin/SRR15859344.bam
#bam/*.bam為bam文件地址
#expr_matrix.txt文件就是下游分析的表達(dá)鉅陣文件
經(jīng)過上面的步驟我們就可以拿到RNA-seq的表達(dá)矩陣了凌外,后面就可以進(jìn)行下游的R語言分析,如差異基因分析涛浙,GO和KEGG分析等康辑。
RNA-seq入門
RNA-seq入門(一)數(shù)據(jù)下載和格式轉(zhuǎn)換
RNA-seq入門(二)質(zhì)控及fastqc報告解讀
RNA-seq入門(三)比對、輸出表達(dá)矩陣