转录组—上游分析GSE232120_微生物物种__NCBI下载参考基因组

原创

sheldor没耳朵

发布于 2024-09-07 10:52:09

860

发布于 2024-09-07 10:52:09

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

转录组—上游分析GSE232120_微生物物种__NCBI下载参考基因组

本次实战的是Hosts Manipulate Lifestyle Switch and Pathogenicity Heterogeneity of Opportunistic Pathogens in the Single-cell Resolution这篇文章。这不是常见的人类或小鼠物种，而是微生物，在上下游处理的过程中与人类或小鼠有所差异。其次我习惯于在ensembl上下载基因组及注释文件，所以首先在这里面查找，但是未找到Serratia marcescens这个物种。本次从NCBI中下载参考基因组及其注释文件。

1 下载raw data

1.1 使用prefetch下载SRR数据

首先将SRR_Acc_List .txt下载至服务器中

## 记得激活环境
cat SRR_Acc_List.txt | while read id
do
echo prefetch ${id} -O ./
done > prefetch.command

less prefetch.command
nohup bash prefetch.command 1>prefetch.log 2>&1 &

1.2 SRR数据转化为fastq数据

cat sra_data/SRR_Acc_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 4 -f fq_data/${id}_1.fastq"
echo "pigz -p 4 -f fq_data/${id}_2.fastq"
done > sra2fq.sh

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

2 fastqc质控

cd /home/data/t110558/project/GSE232120/data/raw_data/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 &

3 trim_galore过滤

cd /home/data/t110558/project/GSE232120/data/raw_data/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

4 Hisat2比对

4.1 微生物Serratia marcescens 参考基因组下载

因为习惯在ensembl上下载基因组及注释文件，所以首先在这里面查找，但是未找到Serratia marcescens这个物种。

于是还是去NCBI里寻找，发现竟然有2712个基因组，选择最新的一个，进入ftp页面即可以看到参考基因组及其注释的下载链接https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/030/291/735/GCF_030291735.1_ASM3029173v1/。

cd /home/data/t110558/database
mkdir GCF_030291735.1
cd GCF_030291735.1
#下载基因组
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/030/291/735/GCF_030291735.1_ASM3029173v1/GCF_030291735.1_ASM3029173v1_genomic.fna.gz >dna.log 
#下载注释文件
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/030/291/735/GCF_030291735.1_ASM3029173v1/GCF_030291735.1_ASM3029173v1_genomic.gff.gz >gff.log

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/030/291/735/GCF_030291735.1_ASM3029173v1/GCF_030291735.1_ASM3029173v1_genomic.gtf.gz >gtf.log
#基因组解压
gunzip GCF_030291735.1_ASM3029173v1_genomic.fna.gz

4.2 参考基因组建立索引

cd $HOME/database/GCF_030291735.1
vim Hisat2Index.sh

mkdir Hisat2Index
fa=GCF_030291735.1_ASM3029173v1_genomic.fna
fa_baseName=GCF_030291735.1.dna
hisat2-build -p 50 -f ${fa} Hisat2Index/${fa_baseName}

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

果然5M的基因组建立索引太快了

4.3 比对参考基因组

时间很长，约12h。

注：pkill -u终止当前用户所有任务

注意在循环中，一般不使用nohup提交命令，不然容易把服务器干爆！因为nohup是把命令放在后台运行。所以每一次循环都会把命令放在后台，导致后台运行命令过多。如果我有5个样本，每个样本用20个线程比对，那一次性就会用100个cpu。或者选择较小的线程，用nohup提交

cd /home/data/t110558/project/GSE232120/data/cleandata

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

4.4 sam转bam

# sam文件转bam
ls *.sam|while read id ;do (nohup samtools sort -O bam -@ 8 -o $(basename ${id} ".sam").bam   ${id}   1>sam2bam.log 2>&1 & );done
# 这个过程会输出大量中间文件
rm *.sam 
# 为bam文件建立索引
ls *.bam |xargs -i samtools index {}

5 FeatureCounts定量

使用featureCounts对bam文件进行定量。

cd ~/project/GSE232120/Expression/featureCounts

## 定义输入输出文件夹
gtf=/home/data/t110558/database/GCF_030291735.1/GCF_030291735.1_ASM3029173v1_genomic.gtf
inputdir=/home/data/t110558/project/GSE232120/Mapping/hisat2

# featureCounts对bam文件进行计数
nohup featureCounts -T 20 -p --countReadPairs -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.bam  1>featureCounts.log 2>&1 &

# 得到表达矩阵
# 处理表头
#/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/GSE232120/Mapping/hisat2/##g' |sed 's#.hisat.bam##g' >raw_counts.txt
#查看表达矩阵
head -n 100 raw_counts.txt|column -t | less -S