生信马拉松 Day18 转录组RNA-seq-3
原创
转录组上游的内容终于上完了，今天的内容太抽象了，每一步处理的内容都不是很好理解，现在上完课也还是摸不着头脑，最大的收获似乎是多按tab键？
可能还是需要多实践来理解
今天的内容主要是比对和定量，用到的包是hisat2、samtools、featurecounts
内容一：hisat2和samtools
比对：质控后数据比对到参考基因组
需要的文件：基因组文件 fasta + 参考数据库注释文件 gtf/gff
可以在相应的官网上下载然后上传到服务器
在Ensembl数据库上下载时需注意：一般下载primary assembly + 注意fasta和gtf/gff文件的版本一致
或使用代码下载
wget -c ftp://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
#这是104版本的，其他版本需要去官网找对应链接，然后修改
wget -c ftp://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz
#这是对应的注释文件，注意版本一致性
# 检查大文件完整性
gzip -t *.gz
# 或者可使用md5码检查Ensembl基因组数据库基因特点：ENSspecies prefix[11位unique digit number】
E exon
G gene
T transcript
例如ENSMUST####，就是Ensembl数据库基因的样式
如果没说物种就是指人类
hisat2用法
hisat2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38.dna
# 这是生成索引的代码，时间较长，最好写成sh脚本提交后台
# 赋值索引前缀
index=/teach/database/GRCh38.104/Hisat2Index/GRCh38.dna
# 单个样本比对
hisat2 -p 2 -x ${index} \
-1 ./trim_galore/SRR1039510_1_val_1.fq.gz \
-2 ./trim_galore/SRR1039510_2_val_2.fq.gz \
-S ./hisat2/SRR1039510.Hisat_aln.sam
# sam 转 bam，并且进行 sort
samtools sort -@ 2 ./hisat2/SRR1039510.Hisat_aln.sam -o ./hisat2/SRR1039510.Hisat_aln.sorted.bam
# 通常来说，不需要生成 sam 文件的，上面几句代码可以通过管道符 | 合并为一句，最后需要有占位符 -
hisat2 -p 2 -x ${index} \
-1 ./trim_galore/SRR1039511_1_val_1.fq.gz \
-2 ./trim_galore/SRR1039511_2_val_2.fq.gz \
| samtools sort -@ 2 -o ./hisat2/SRR1039511.Hisat_aln.sorted.bam -
# 对bam建索引，这是一句代码
samtools index ./hisat2/SRR1039510.Hisat_aln.sorted.bam ./hisat2/SRR1039510.Hisat_aln.sorted.bam.baisam保留了序列比对文件，是比对后的文件
数据变换的过程是：sra-fastq-sam/bam
bam是sam的二进制文件
samtools flags 99
显示0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1
0x63是一个数字取和，类似linux里给权限的取值求和，0x63=0x1+0x2+0x20+0x40，分别对应PAIRED,PROPER_PAIR,MREVERSE,READ1，此时也叫99
这个内容有对应的表
samtools view 主要是看bam文件，因为sam可以用zless等方式都可以看，而bam是压缩文件
内容二：featureCounts
featurecount 表达定量
得到的结果是原始表达矩阵raw counts然后再处理才得到clean data
## 定义输入输出文件
gtf=/teach/database/GRCh38.104/Homo_sapiens.GRCh38.104.chr.gtf.gz
featureCounts -a $gtf -o ./featureCounts/all.count.txt -p -T 6 -t exon -g gene_id ./hisat2/*.sorted.bam
# 得到表达矩阵
# 处理表头，要换成自己的路径
grep -v '#' ./featureCounts/all.count.txt |cut -f 1,7- | sed "s@./hisat2/@@g" | sed 's#.Hisat_aln.sorted.bam##g' > ./featureCounts/raw_counts.txt
sed s///
sed s##
sed s%%%
# 列对齐显示
head ./featureCounts/raw_counts.txt | column -t生信技能树，生信马拉松，火龙果老师
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。