有参转录组的比对

一、数据准备

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 “$3” -S “$2".sam --rna-strandness R 1>"$2".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 "$2".bam "$2".sam &"}' samples.txt的路径 >samtools.sh
nohop sh samtools.sh &

:wqcd

©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容