有参转录组的比对

一、数据准备

1、测序数据

（1 ）公司获得
（2） NCBI下载/ENA下载

2、参考序列准备

（1）基因组序列 下载地址：Ensembl、EnsemblGenomes、NCBI
下载后解压整合到一起构成genome.fasta文件
（2）基因注释文件 下载地址：Ensembl、EnsemblGenomes、NCBI 如果下载不到gtf,而是gff文件，可以用gffread命令转换

gffread -T -o 输出文件名 输入文件路径

3 测序数据的样本名及测序数据的绝对路径

如果是双末端测序，为4列，第一列为样本分组名 第二列为样本名 第三列为测序数据的绝对路径

cat sample.txt
ck_b001   ck_b001_s1  /home/luly/workspace/data/ck_b001_s1.fq.gz
ck_b001   ck_b001_s2  /home/luly/workspace/data/ck_b001_s2.fq.gz
ck_b001   ck_b001_s3  /home/luly/workspace/data/ck_b001_s3.fq.gz

这个表可以用来批量生成脚本，所以要非常仔细

二、比对到参考基因组

1、安装软件

hisat2
samtools

2、为参考序列构建索引

hisat2-built ../ref/genome.fasta ../ref/genome 1>histat2-built.log 2>&1

hisat2-built 基因组文件（genome.fasta）路径 索引名称路径 1>histat2-built.log 2>&1

3、批量比对

awk '{print "hisat2 -- new-summary -p 10 -x 索引名称路径 -U “” -S “".sam --rna-strandness R 1>"".log 2>&1 &"}' samples.txt的路径 >hisat.sh
#-p 线程数
#-U 单端测序文件，fastq格式。双端的为-1 -2
#-S 输出文件，sam格式
#--rna-strandness  普通文库，单端为R，双端为RF
#1>"&2".log 日志文件
#2>&1 错误输出也定向到log文件
#最后一个&符号是让它并行
nohop sh hisat.sh &
#后台运行

4、比对结果压缩(批量）

由于sam文件太大。需要将sam文件压缩为bam文件并排序

awk '{print "samtools sort -o "".bam "".sam &"}' samples.txt的路径 >samtools.sh
nohop sh samtools.sh &

:wqcd