1 修改后缀

如果是数据库公开数据，从此步骤开始：

cd /work/home/acwnw4bl7y/HT22/HT22
ls *lite.1|while read id;
do
mv $id ${id%%.*}.sra;
done

2 数据解压缩

cd /work/home/acwnw4bl7y/HT22/HT22
export PATH=/work/home/acwnw4bl7y/sratoolkit.3.1.1-centos_linux64/bin:$PATH 
ls *.sra|while read id;
do
fastq-dump --gzip --split-3 $id;
done

cd /work/home/acwnw4bl7y/DRG1
source /work/home/acwnw4bl7y/apprepo/sra_tools/3.1.1-gnu7.3.1/scripts/env.sh
ls *.sra|while read id;
do
fastq-dump --gzip --split-3 $id;
done

3 修剪接头序列

如果是自己测序产生的数据，从此步骤开始：

cd /work/home/acwnw4bl7y/HT22/HT22
source ~/miniconda3/bin/activate
conda activate pyht
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

cd /work/home/acwnw4bl7y/DRG1
source /work/home/acwnw4bl7y/apprepo/fastp/0.19.10-gcc731/scripts/env.sh
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 比对基因组

cd /work/home/acwnw4bl7y/HT22/HT22
source ~/miniconda3/bin/activate
conda activate pyht
ls *1.clean.fq.gz|while read id;
do
hisat2 -t -p 16 -x /work/home/acwnw4bl7y/GRCm39/GRCm39 \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.hisat2.log \
|samtools sort -@ 16 -o ${id%_*}_ht2p.bam
done

cd /work/home/acwnw4bl7y/DRG1
source /work/home/acwnw4bl7y/apprepo/hisat2/2.2.1-none/scripts/env.sh
source /work/home/acwnw4bl7y/apprepo/samtools/1.19.2-none/scripts/env.sh
ls *1.clean.fq.gz|while read id;
do
hisat2 -t -p 16 -x /work/home/acwnw4bl7y/GRCm39/GRCm39 \
-1 $id -2 ${id%_*}_2.clean.fq.gz \
2>${id%%_*}.hisat2.log \
|samtools sort -@ 16 -o ${id%_*}_ht2p.bam
done

4 Count计数

cd /work/home/acwnw4bl7y/HT22/HT22
source ~/miniconda3/bin/activate
conda activate pyht
mkdir -p htseq
ls *.bam |while read id;
do 
htseq-count -f bam -s no -i gene_id $id \
/work/home/acwnw4bl7y/Mus_musculus.GRCm39.106.gtf \
1>./htseq/${id%%_*}.txt 2>./htseq/${id%%_*}.HTseq.log;
done

cd /work/home/acwnw4bl7y/DRG1
module load /work/home/acwnw4bl7y/apprepo/htseq/v2.0.1-none/scripts/htseq-v2.0.1-none
mkdir -p htseq
ls *.bam |while read id;
do 
htseq-count -f bam -s no -i gene_id $id \
/work/home/acwnw4bl7y/Mus_musculus.GRCm39.106.gtf \
1>./htseq/${id%%_*}.txt 2>./htseq/${id%%_*}.HTseq.log;
done