转录组上游分析—使用iseq下载原始数据、小鼠基因组、单端测序数据处理

原创

sheldor没耳朵

发布于 2024-08-21 11:16:22

1940

发布于 2024-08-21 11:16:22

文章被收录于专栏：数据挖掘

转录组上游分析—使用iseq下载原始数据、小鼠基因组、单端测序数据处理

1 下载rawdata_使用iseq替代prefetch进行数据下载

进行数据集GSE105789上游分析的时候，总共才四个数据集，使用prefetch下载的时候，不知道网络抽了什么风，速度一直都很慢。下了10个小时才下了三分之一。!

经过大佬指点，尝试使用iseq替代prefetch进行数据下载

参考链接：https://mp.weixin.qq.com/s/AVqv07swFvjl6OCnLwwLPA

iseq可以直接下载GSE编号，不用再使用prefetch下载list中的SRA号，而且不需要再通过fastq-dump进行从SRA到fastq数据的转换。

首先我在当前conda环境中直接安装iseq是有冲突的，因此需要重新建立个conda环境：

正式使用

conda create -n iseq
conda activate iseq
conda install iseq
#直接传输GSE编号即可，不用再使用prefetch下载list中的SRA号
iseq -i GSE105789 -g

开始还是很快的，可以根据success.log中的内容确定哪些文件是已经成功下载的，

如果终止了重新下载需要把未下载完全的文件删除掉，如图中如果需要重新下载，需要删除掉SRR6203123,

同时可以看到下载下来的直接就是fastq.gz的文件，不必再进行sra to fastq的转换。

如图，重新下载会跳过已经下载好的。

后面的速度有所下降，可能是网络波动的问题，但是总体还是比prefetch快的。

这里还是采用prefetch下载(因为之前已经下载了90%)，然后通过fasterq-dump进行转换

单端数据

cat sra_data/SRR_download_List.txt | while read id
do
echo "fasterq-dump -e 32 --split-files -O fq_data/ --outfile ${id}.fastq sra_data/${id}/${id}.sra"    
# 这一步使用fastq-dump或fasterq-dump都可以
echo "pigz -p 16 -f fq_data/${id}.fastq"
done > sra2fq.sh

less sra2fq.sh
nohup bash sra2fq.sh 1>sra2fq.log 2>&1 &

2 fastqc质控\trim_galore过滤

cd /home/data/t110558/project/GSE105789/data/rawdata/qc
nohup fastqc -t 20 -o ./ ../fq_data/SRR*.fastq.gz 1>qc.log 2>&1 &
#数据整合
multiqc *.zip -o ./ -n qc_fastqc  1>./multiqc.log 2>&1 &

#样本数量少就可以直接提交到后台运行
cd /home/data/t110558/project/GSE105789/data/rawdata/fq_data/
ls *gz|cut -d"." -f 1|sort -u |while read id;do
nohup trim_galore -j 20 -q 25 --phred33 --length 36 --stringency 3 --paired -o ../../cleandata   ${id}*.gz 1>../../raw2clean.log 2>&1 & 
done

#修改代码
cd /home/data/t110558/project/GSE105789/data/rawdata/fq_data/
ls *gz|cut -d"." -f 1|sort -u |while read id;do
nohup trim_galore -j 8 -q 25 --phred33 --length 36 --stringency 3  -o ../../cleandata   ${id}*.gz 1>../../raw2clean.log 2>&1 & 
done

3 hisat2比对

注意要使用小鼠的基因组

cd /home/data/t110558/project/GSE105789/data/cleandata

indexPrefix=/home/data/t110558/database/GRCh38.111/Hisat2Index/GRCh38.dna
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup  hisat2 -p 20 -x  $indexPrefix -U ${id}*_trimmed.fq.gz  -S ../../Mapping/hisat2/${id}.hisat.sam 1>../../Mapping/hisat2/hisat2.log 2>&1 & 
done

参考基因组选择小鼠

3.1 小鼠参考基因组下载

# 下载基因组序列
nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz  >dna.log &

# 下载基因组注释文件
nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/gtf/mus_musculus/Mus_musculus.GRCm39.112.gtf.gz >gtf.log &

nohup axel -n 20 https://ftp.ensembl.org/pub/release-112/gff3/mus_musculus/Mus_musculus.GRCm39.112.gff3.gz >gff.log&
#解压对应文件
nohup gunzip Mus_musculus.GRCm39.112.gtf.gz Mus_musculus.GRCm39.112.gff3.gz Mus_musculus.GRCm39.dna.primary_assembly.fa.gz >unzip.log &

3.2 参考基因组建立索引

cd /home/data/t110558/database/GRCm39.112
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
vim Hisat2Index.sh

mkdir Hisat2Index
fa=Mus_musculus.GRCm39.dna.primary_assembly.fa
fa_baseName=GRCm39.dna
hisat2-build -p 50 -f ${fa} Hisat2Index/${fa_baseName}

# 运行
nohup bash Hisat2Index.sh >Hisat2Index.sh.log &

3.3 基因组比对

cd /home/data/t110558/project/GSE105789/data/cleandata

indexPrefix=/home/data/t110558/database/GRCm39.112/Hisat2Index/GRCm39.dna
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup  hisat2 -p 20 -x  $indexPrefix -U ${id}*_trimmed.fq.gz  -S ../../Mapping/hisat2/${id}.hisat.sam 1>../../Mapping/hisat2/hisat2.log 2>&1 & 
done

#sam转bam
ls *.sam|while read id ;do (nohup samtools sort -O bam -@ 20 -o $(basename ${id} ".sam").bam   ${id}   1>sam2bam.log 2>&1 & );done
#为bam文件建立索引
ls *.bam |xargs -i samtools index {}
rm *.sam

4 FeatureCounts定量

注意是单端数据

cd /home/data/t110558/project/GSE105789/Expression/featureCounts

## 定义输入输出文件夹
gtf=/home/data/t110558/database/GRCm39.112/Mus_musculus.GRCm39.112.gtf
inputdir=/home/data/t110558/project/GSE105789/Mapping/hisat2

# featureCounts对bam文件进行计数
nohup featureCounts -T 20  -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.bam  1>featureCounts.log 2>&1 &
# 对定量结果质控
multiqc all.id.txt.summary

# 得到表达矩阵
# 处理表头
#/home/data/t110558/project/GSE105789/Mapping/hisat2/SRR6203120.hisat.bam
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#/home/data/t110558/project/GSE105789/Mapping/hisat2/##g' |sed 's#.hisat.bam##g' >raw_counts.txt
#查看表达矩阵
head -n 100 raw_counts.txt|column -t