cd /pnas/fangxd_group/renyx/macaque/01rawdata for ((i=230;i<=237;i++)) ;do /software/biosoft/software/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump --split-3 --gzip SRR4042$i.sra;done
质控
mkdir QC echo "fastqc started at $(date)" #输出运行时间 fastqc -o QC *gz multiqc *fastqc.zip --ignore *.html echo "fastqc finished at $(date)"
质控检测采用fastqc和multiqc联合使用, multiqc有利于多个样本同时展示质量检测结果。FastQC的检测报告左侧会给出测序质量的一个summary,通常并不是所有的检测项都会是绿色通过,会有一些警告和fail的项目。主要可以从Per base sequence quality 看一下测序碱基质量,Per sequence GC content 看一下GC含量,如果实际的GC含量(红线)出现双峰,且导致后期的序列比对很低时,需要引起注意,有可能是存在样品污染;再者就是看一下Adapter是否去除及去除的是否干净。这里的Adapter虽然显示通过,但是去除的并不干净,所以后续还会进一步去除Adapter。
去除接头和低质量值
Trimmomatic 支持多线程,处理数据速度快,主要用来去除 Illumina 平台的 fastq 序列中的接头,并根据碱基质量值对 fastq 进行修剪。软件有两种过滤模式,分别对应 SE 和 PE 测序数据,同时支持 gzip 和 bzip2 压缩文件。另外也支持 phred-33 和 phred-64 格式互相转化。
#下载subread wget https:///projects/subread/files/subread-1.5.3/subread-1.5.3-Linux-x86_64.tar.gz tar zxvf subread-1.5.3-Linux-x86_64.tar.gz
常用参数说明:
-p 针对paired-end数据 T 多线程数 t 默认将exon作为一个feature g 默认是gene_id,从注释文件中提取Meta-features信息用于read count a 基因组注释文件,默认是gtf o 输出文件
(PS:存储不够了,后续选择两组数据做流程分析。)
echo"hisat2 started at $(date)" #make index cd /pnas/fangxd_group/renyx/macaque/00ref hisat2-build -p 8 Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa maca_index
#alignment for i in {230..231} do hisat2 -t -x 00ref/maca_index -1 02clean_data/SRR4042$i\_paired_clean_R1.fastq.gz \ -2 02clean_data/SRR4042$i\_paired_clean_R2.fastq.gz \ -S 03align_out/SRR4042$i\.sam done
#covert to bam, and sort bam for i in {230..231} do samtools view -Sb SRR4042$i\.sam > SRR4042$i\.bam samtools sort SRR4042$i\.bam -o SRR4042$i\_sorted.bam samtools index SRR4042$i\_sorted.bam rm *sam done echo"hisat2 finished at $(date)"
#featureCounts定量 cd /pnas/fangxd_group/renyx/macaque/ featureCounts=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/featureCounts $featureCounts -p -T 6 -t exon -g gene_id -a 00ref/Macaca_mulatta.Mmul_8.0.1.91.gtf \ -o 05featurecount_out/counts.txt 03align_out/hisat2_mapping/*sorted.bam
#make subread index echo"subread started at $(date)" buildindex=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subread-buildindex cd /pnas/fangxd_group/renyx/macaque/00ref $buildindex -o maca Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa
#alignment subjunc=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subjunc subjunc_maca_index=/pnas/fangxd_group/renyx/macaque/00ref/maca cd /pnas/fangxd_group/renyx/macaque/ for i in {230..231} do $subjunc -T 5 -i $subjunc_maca_index -r 02clean_data/SRR4042$i\_paired_clean_R1.fastq.gz -R 02clean_data/SRR4042$i\_paired_clean_R2.fastq.gz -o 03align_out/subread_mapping/SRR4042$ done