hisat2 -t -x ~/rna_seq_analysis/reference/index/hg19/genome -1 ~/rna_seq_analysis/fastq/SRR5894154_1.fastq.gz -2 ~/rna_seq_analysis/fastq/SRR5894154_2.fastq.gz -S ~/rna_seq_analysis/aligned/SRR5894154.sam
之前出问题好像是因为少空了空格。
比对结果
http://blog.sciencenet.cn/blog-3334560-1078097.html
得到的sam格式文件40多G。。。。
我哭了。。。
下一步把sam文件转化为bam文件，用samtools
SAM格式是目前用来存放大量核酸比对结果信息的通用格式，也是人类能够“直接”阅读的格式类型，而BAM和CRAM是为了方便传输，降低存储压力将SAM进行压缩得到的格式形式。 注，BAM格式必须要建立索引才能快速读取指定位置的信息。
# 1.3版本前
samtools view -bS bwa.sam > bwa.bam
samtools sort bwa.bam > bwa_sorted.bam
samtools index bwa_sorted.bam
# 1.3版本后
samtools sort bwa.sam > bwa_sorted.bam
samtools index bwa_sorted.bam
1. 格式转换
2. 排序
3. 索引
大于号：将一条命令执行结果（标准输出，或者错误输出，本来都要打印到屏幕上面的）重定向其它输出设备（文件，打开文件操作符，或打印机等等）
1. samtools view -S SRR5894154.sam -b > SRR5894154.bam
bam文件不到8g，于是赶紧把sam文件删了～
2. samtools sort SRR5894154.bam > SRR5894154_sorted.bam
3. samtools index SRR5894154_sorted.bam
//虽然我不知道第三步有什么用。。。
如果你要比较同一个样本(within-sample)不同基因之间的表达情况，你就需要考虑到转录本长度，因为转录本越长，那么检测的片段也会更多，直接比较等于让小孩和大人进行赛跑。如果你是比较不同样本（across sample）同一个基因的表达情况，虽然不必在意转录本长度，但是你要考虑到测序深度（sequence depth)，毕竟测序深度越高，检测到的概率越大。除了这两个因素外，你还需要考虑GC%所导致的偏差，以及测序仪器的系统偏差。目前对read count标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之间的差异与换算方法已经有文章进行整理和吐槽了。
在转录本水平上，一般常用工具为Cufflinks和它的继任者StringTie， eXpress。这些软件要处理的难题就时转录本亚型（isoforms）之间通常是有重叠的，当二代测序读长低于转录本长度时，如何进行区分？这些工具大多采用的都是expectation maximization（EM）。好在我们有三代测序。上述软件都是alignment-based，目前许多alignment-free软件，如kallisto, silfish, salmon，能够省去比对这一步，直接得到read count，在运行效率上更高。不过最近一篇文献[1]指出这类方法在估计丰度时存在样本特异性和读长偏差。
-f bam/sam： 指定输入文件格式，默认SAM
-r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序，默认是name
-s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的，所以需要判断到底来自于哪条链。如果选择了no， 那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上，第二个read则在另一条链上。
-a 最低质量， 剔除低于阈值的read
-m 模式 union（默认）, intersection-strict and intersection-nonempty。一般而言就用默认的，作者也是这样认为的。
-i id attribute: 在GTF文件的最后一栏里，会有这个基因的多个命名方式（如下）， RNA-Seq数据分析常用的是gene_id， 当然你可以写一个脚本替换成其他命名方式。
gene_id "ENSG00000223972.5_2"; transcript_id "ENST00000456328.2_1"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1_1"; level 2;.
htseq-count -s no -r pos -f bam ~/rna_seq_analysis/aligned/SRR5894154_sorted.bam ~/rna_seq_analysis/human_genome/gencode.v28lift37.annotation.sorted.gff3 > ~/rna_seq_analysis/aligned/SRR5894154.count
网友评论