这个部分主要是序列比对和reads读取,别人的帖子写的很全面,但有点复杂,需要兼顾来看。
//www.greatytc.com/p/ff585b72f04e
http://www.biotrainee.com/home.php?mod=space&uid=424&do=thread&view=me&from=space
总体来说,就是首先,在HIST网站下载人类的序列(我理解就是index),可以用终端下,也可以在HIST网站下,在网站下会快很多,然后放到终端解压。然后就是之前转换好的fastq格式(有1和2的那个),经过index比对,形成sam文件,详细过程看下面命令,每个编号转换为sam大概35分钟左右(根据大小和计算机水平),这个sam真心大,随便一个都十几个G。转换完sam之后,需要把它转换为bam格式,然后序列重整之后再转回sam.count格式。然后就可以reads了,经过reads之后就会变成count了,就可以导入R里面做差异基因等进一步分析了。
比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。 直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。 接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看! 顺便对bam文件进行简单QC,参考直播我的基因组系列。
下载 index
axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz #解压
ps.可以再HIST下,对应物种,老鼠就下老鼠的。。。
正式比对
命令行目录:Desktop/zhuanlu_files
hg19目录:Desktop/zhuanlu_files/hg19
RNA-Seq Data: Desktop/zhuanlu_files/RNA-Seq
注释:
-t 记录时间
-x hg19(index)文件路径
-1 -2 测序的两个fastq文件(单端测序用-U)
-S 比对结果输出路径
mkdir -p RNA-Seq/aligned
for i in`seq 56 58`
do
hisat2 -t -x hg19/genome -1 RNA-Seq/SRR35899${i}.sra_1.fastq.gz -2 RNA-Seq/SRR35899${i}.sra_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam
done
上步就是把fastq格式转换为sam格式,并列序列比对
比对解读
总共比对了28856780条,7.4%一次都没有比对到,77.18%比对到一次,15.41%大于一次
没有匹配上的7.4%,不按照顺序再匹配,有4.3%匹配了
之后把剩余部分用单端数据进行比对,64.28%没有比对上,26.16比对了一次,9.56%大于一次
总共比对到95.45%
SAMtools
SAMtools 详解:
//www.greatytc.com/p/15f3499a6469
SAM 数据格式是目前高通量测序中存放比对数据的标准格式,但是SAM 文件都很大,非常占空间,所以需要转到bam文件,而用的就是SAMtools软件。
注释
-S 输入sam文件
-b 指定输出的文件为bam
排序默认是根据染色体的位置
for i in`seq 59 62`
do
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam
samtools view -h SRR35899${i}_nsorted.bam > SRR35899${i}_count.sam
done
上步就是先将sam 文件转成bam文件,bam 文件用 read name (-n)排序,再转回到sam文件,把这些转换好的文件放上层文件夹也行,新建文件夹也行,直接放上层文件夹会方便一些。
mkdir -p RNA-Seq/matrix/
htseq-count ~/SRR/SRR3589918_count.sam ~/gencode.v27lift37.annotation.gtf > 18.count
Genecode/gencode.v27lift37.annotation.gtf 这个是之前,在Genecode网站上下好的注释,老鼠的注释可以看看原文。把gtf这个文件放入Genecode这个文件夹,下面是原文的命令代码。
htseq-count RNA-Seq/aligned/SRR3589959_count.sam mouse/gencode.vM10.annotation.gtf
htseq-count RNA-Seq/aligned/SRR3589960_count.sam mouse/gencode.vM10.annotation.gtf
htseq-count RNA-Seq/aligned/SRR3589961_count.sam mouse/gencode.vM10.annotation.gtf
htseq-count RNA-Seq/aligned/SRR3589962_count.sam mouse/gencode.vM10.annotation.gtf
之后就可以得到count的文件,进行下一步处理。