1 软件安装
https://www.jianshu.com/p/eb89ab4af035
linux平台下需要安装的软件:fastqc,fastp,hisat2,samtools,htseq
2 线虫基因组和基因组注释文件的下载
http://asia.ensembl.org/index.html
http://asia.ensembl.org/Caenorhabditis_elegans/Info/Index
基因组:
wget -c ftp://ftp.ensembl.org/pub/release-101/fasta/caenorhabditis_elegans/dna/Caenorhabditis_elegans.WBcel235.dna.toplevel.fa.gz
关于nonchromosomal,toplevel,rm.nonchromosomal,rm.toplevel,sm.nonchromosomal,sm.toplevel的区别参考:
关于人类参考基因组的一些认识
下载基因组注释文件:
wget -c ftp://ftp.ensembl.org/pub/release-101/gtf/caenorhabditis_elegans/Caenorhabditis_elegans.WBcel235.101.gtf.gz
构建比对索引文件
将基因组文件解压为fasta文件(Caenorhabditis_elegans.WBcel235.dna.toplevel.fa)
hisat2-build -p 2 Caenorhabditis_elegans.WBcel235.dna.toplevel.fa Caenorhabditis_elegans
索引文件也可以再Hisat2官网直接下载:这样就不用下载基因组序列文件了。
wget -c https://genome-idx.s3.amazonaws.com/hisat/wbcel235.tar.gz
3 下机数据的质控和过滤
3.1 质控
mkdir -p fastqc
ls *.fq.gz|while read id;
do
fastqc -t 2 $id -o ./fastqc;
done
3.2 过滤
mkdir -p fastp
ls *1.fq.gz|while read id;
do
fastp -5 20 -i ${id%_*}_1.fq.gz -I ${id%_*}_2.fq.gz \
-o ${id%_*}_1.clean.fq.gz -O ${id%_*}_2.clean.fq.gz \
-j ./fastp/${id%_*}.json -h ./fastp/${id%_*}.html;
done
4 比对以及转换格式
ls *1.clean.fq.gz|while read id;
do
hisat2 -t -p 2 -x /media/luozhixin/0000678400004823/Indexs/Hisat2/Caenorhabditis_elegans/Caenorhabditis_elegans \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.ht2.log \
|samtools sort -@ 2 -o ${id%_*}_ht2p.bam
done
5 生成reads count
ls *.bam |while read id;
do
htseq-count -f bam -s no -i gene_id $id \
/media/luozhixin/0000678400004823/Gtf_gff/Caenorhabditis_elegans/Caenorhabditis_elegans.WBcel235.101.gtf \
1>${id%%.*}.txt 2>${id%%.*}.HTseq.log;
done
6 组成count计数矩阵
利用python脚本,合并所有样本的reads count,也可以在Windows下完成,需要安装windows版本的python
下载注释文件:
Caenorhabditis elegans (ID 41) - Genome - NCBI (nih.gov)
Genome List - Genome - NCBI (nih.gov)
网友评论