1. 下载数据
高通量测序数据检索及下载参考：//www.greatytc.com/p/6798540d2294
示例
prefetch SRR8786111 -O ./
prefetch SRR8786110 -O ./
或者
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/
2. 格式转换
ls *.sra|while read id;
do
fastq-dump --gzip --split-3 $id;
done
3. 质控
mkdir -p fastqc
ls *.fq.gz|while read id;
do
fastqc -t 2 $id -o ./fastqc;
done
t，线程
4. 过滤
4.1 双末端过滤
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
如果是SRA数据库下载的数据：
mkdir -p fastp
ls *1.fastq.gz|while read id;
do
fastp -5 20 -i ${id%_*}_1.fastq.gz -I ${id%_*}_2.fastq.gz \
-o ${id%_*}_1.clean.fq.gz -O ${id%_*}_2.clean.fq.gz \
-j ./fastp/${id%_*}.json -h ./fastp/${id%_*}.html;
done
4.2 单末端过滤
mkdir -p fastp
ls *.fq.gz|while read id;
do
fastp -5 20 -i $id -o ${id%%.*}.clean.fq.gz \
-h ./fastp/${id%%.*}.html -j ./fastp/${id%%.*}.json;
done
-a，接头序列；
-5，5’端测序质量值；
--length_required，最短长度
--length_limit，最长长度
如果是SRA数据库下载的数据：
mkdir -p fastp
ls *.fastq.gz|while read id;
do
fastp -5 20 -3 20 -i $id -o ${id%%.*}.clean.fq.gz \
-h ./fastp/${id%%.*}.html -j ./fastp/${id%%.*}.json;
done
5.比对及SAM格式转换
此处使用了管道命令
5.1 下游为7，8
ls *1.clean.fq.gz|while read id;
do
hisat2 --dta -t -p 2 -x /home/luozhixin/grch38_tran/genome_tran \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.hisat.log \
|samtools sort -@ 2 -o ${id%_*}_ht2p.bam
done
p,线程，
x,参考基因组
samtools sort：这一步主要是把空间占用较大的sam文件转换为空间占用较小的二进制文件，即bam文件，并将其按比对到的染色体位置或read的ID排序。
-n按reader ID sort，默认按染色体位置sort。
注意：以染色体位置sort要比reader ID sort节省空间
5.2 下游为8
双末端
ls *1.clean.fq.gz|while read id;
do
hisat2 -t -p 2 -x /home/luozhixin/hg19/genome \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.hisat.log \
|samtools sort -@ 2 -o ${id%_*}_ht2p.bam
done
单末端
ls *.clean.fq.gz|while read id;
do
hisat2 -t -p 2 -x /media/luozhixin/本地磁盘/bioinfomatics/INDEX/hisat2/human/grch38/genome \
-U $id \
2>${id%%.*}.hisat2.log \
|samtools sort -@ 2 -o ${id%%.*}_ht2p.bam
done
p线程，x参考基因组
7. Stringtie
7.1 初始转录本组装
mkdir -p htseq
ls *.bam |while read id;
do
stringtie $id -G /home/luozhixin/gencode.v28lift37.annotation.gtf -o ${id%_*}.gtf -p 2 #G，gtf p线程
echo -e "GTF/${id%.*}.gtf" >> GTF/mergelist.txt
done
7.2 merge转录本
ls *.gtf |while read id;
do
stringtie --merge -p 2 -G /home/luozhixin/gencode.v28lift37.annotation.gtf -o stringtie_merge.gtf GTF/mergelist.txt #G，gtf p线程
done
7.3 重新组装转录本
ls *.bam |while read id;
do
stringtie $id -G stringtie_merge.gtf -o ballgown/${id%_*}/${id%_*}.gtf -A ${id%_*}.tab -B -e -p 2 #G，gtf p线程
done
7.4 转录本计数
prepDE.py -g gene_count_matrix.csv -t transcript_count_matrix.csv
8. HTseq
根据情况选择以下任意一种计数方法
8.1 非链特异性建库
8.1.1 bam文件以pos进行sort情况下：
ls *.bam |while read id;
do
htseq-count -f bam -r pos -s no -i gene_id $id \
/home/luozhixin/gencode.v28lift37.annotation.gtf 1>${id%%.*}.txt 2>${id%%.*}.HTseq.log;
done
8.1.2 bam文件以name进行sort情况下：
ls *.bam |while read id;
do
htseq-count -f bam -s no -i gene_id $id \
/home/luozhixin/gencode.v28lift37.annotation.gtf 1>${id%%.*}.txt 2>${id%%.*}.HTseq.log;
done
8.2 链特异性建库
mkdir -p htseq
ls *.bam |while read id;
do
htseq-count -f bam -s reverse -t gene -i gene_id $id \
/home/luozhixin/gencode.v28lift37.annotation.gtf \
1>./htseq/${id%%.*}.txt 2>./htseq/${id%%.*}.HTseq.log;
done
f，指定输入bam或sam格式；
s，是否为链特异性建库；
-i，以什么为单位计数；gene_id或gene_name，一般不使用transcript_id或transcript_name
1，指定输出文件名；
2，指定日志文件名
