转录组—上游分析_如何拿到count矩阵

2024-08-12 18:56:06 浏览数 (2)

转录组—上游分析_如何拿到count矩阵

本文档记录GSE149638数据集中下载SRR11652578和SRR11652615原始数据

将其从SRR格式转换为fastq格式，通过fastqc质控，trim_galore过滤，Hisat2比对至基因组，FeatureCounts定量，最终得到基因表达的count矩阵的过程。

1 文件目录

文件目录参考如下，每一个项目就在project中建立相应的文件夹

代码语言：bash复制
## 示例如下：
├── database # 数据库存放目录，包括参考基因组，注释文件，公共数据库等
├── project  # 项目分析目录
    └── Human-16-Asthma-Trans #具体项目
        ├── data # 数据存放目录
        │   ├── cleandata # 过滤后的数据
           	│	├── trim_galore # trim_galore过滤
		   	│	└── fastp	    # fastp过滤
        │   └── rawdata # 原始数据
        ├──  Mapping # 比对目录
        │   ├── Hisat2 # Hisat比对
        │   └── Subjunc # subjunc比对
        └── Expression # 定量
            ├── featureCounts # featureCounts
            └── Salmon # salmon定量
        └── Diff_analysis # 差异分析

2 下载raw data

下载原始数据GSE149638

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA629498&o=acc_s:a

从96个数据集中下载最小的两个数据集作为测试

注：

“prefetch fasterq-dump 的组合是从 SRA 访问中提取 FASTQ 文件的最快方法。预取工具将所有必要的文件下载到您的计算机上。如果下载不成功，可以多次调用 prefetch - 工具。它不会每次都从头开始;相反，它将从上次调用失败的位置开始。https://github.com/ncbi/sra-tools/wiki/08.-prefetch-and-fasterq-dump”

2.1 安装SRA Toolkit

https://github.com/ncbi/sra-tools/wiki/02.-Installing-SRA-Toolkit

代码语言：bash复制
wget --output-document sratoolkit.tar.gz https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
#下载速度太慢，从网页下载再上传至服务器
tar -vxzf sratoolkit.current-ubuntu64.tar.gz 
export PATH=$PATH:$PWD/sratoolkit.3.1.1-ubuntu64/bin
which fastq-dump
#/home/data/t110558/sratoolkit.3.1.1-ubuntu64/bin/fastq-dump
#测试工具包是否正常运行
fastq-dump --stdout -X 2 SRR390728

2.2 使用prefetch下载SRR数据

参考链接https://mp.weixin.qq.com/s/M8-ZuGgW2TQf2cKF4Mea3A

SRR_download_List.txt

代码语言：bash复制
## 记得激活环境
cat SRR_download_List.txt | while read id
do
echo prefetch ${id} -O ./
done > prefetch.command

less prefetch.command
nohup bash prefetch.command &

2.3 SRR数据转化为fastq数据

代码语言：bash复制
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}_1.fastq"
echo "pigz -p 16 -f fq_data/${id}_2.fastq"
done > sra2fq.sh

less sra2fq.sh
nohup bash sra2fq.sh &

注：

生成一个名为 sra2fq.sh 的脚本文件，脚本内容是用来从 .sra文件中提取 .fastq 文件，并压缩这些 .fastq 文件。以下是代码的逐行解释：

  1. cat ../SRR_Acc_List.txt | while read id:
    • cat ../SRR_Acc_List.txt` 命令读取SRR_Acc_List.txt 文件的内容，文件中可能存储了多个SRR ID，每一行一个。
    • while read id 表示逐行读取该文件的内容，并将每一行的内容赋值给变量 id，以便在循环中使用。
  2. do:
    • 开始一个循环体，对于 SRR_Acc_List.txt 中的每个 id，循环体中的命令将会被执行。
  3. echo "fasterq-dump -e 32 --split-files -O ./ --outfile ${id}.fastq ${id}.sra"`:
    • 这行代码生成一个 fasterq-dump 命令，用于将 .sra文件转换为 .fastq 文件。
    • -e 32 指定使用 32 个线程来加速处理。
    • --split-files 参数用于将配对的读段（paired-end reads）分为两个不同的 .fastq` 文件（即 _1.fastq 和 _2.fastq）。
    • -O ./ 指定输出目录为当前目录。
    • --outfile ${id}.fastq 指定输出文件的前缀。
    • ${id}.sra 指定输入的 .sra 文件名。
  4. echo "pigz -p 16 -f ./${id}_1.fastq":
    • 生成一个 pigz 压缩命令，用于压缩上一步生成的 ${id}_1.fastq 文件。
    • -p 16 指定使用 16 个线程来加速压缩。
    • -f强制压缩，即使目标文件已经存在。
  5. echo "pigz -p 16 -f ./${id}_2.fastq":
    • 同样生成一个 pigz 压缩命令，用于压缩 ${id}_2.fastq 文件。
  6. done > sra2fq.sh:
    • 循环结束。
    • 将所有生成的命令行输出重定向到 sra2fq.sh 文件中。这样，sra2fq.sh 文件中将包含针对每个 SRR ID 的一系列命令，用于提取 .fastq文件并进行压缩。

3 fastqc质控

代码语言：bash复制
cd /home/data/t110558/project/Human-16-Asthma-Trans/data/rawdata/qc
nohup fastqc -t 20 -o ./ ../fq_data/SRR*.fastq.gz >qc.log &
#数据整合
multiqc *.zip -o ./ -n qc_fastqc  >./multiqc.log

qc_fastqc.html

4 trim_galore过滤

代码语言：bash复制
#样本数量少就可以直接提交到后台运行
cd $HOME/project/Human-16-Asthma-Trans/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 & 
done

#如果数量多,可以用这个脚本控制下并行数量

代码语言：bash复制
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
echo trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ../2.1.clean_fq   ${id}_*.gz 
done   > trim_galore.sh 

for i in {0..1};do ( nohup bash ../submit.sh trim_galore.sh  2 $i 1>log.trim_galore.$i.txt 2>&1 & )  ;done 

# 如果是单端 ： 
ls *gz| sort -u |while read id;do
echo trim_galore -q 25 --phred33 --length 36 --stringency 3   -o ../cleanData   ${id} 
done   > trim_galore.sh

这样就是每次并行10个样品，相当于批处理啦，其中 submit.sh 脚本内容如下所示；

代码语言：bash复制
cat $1 |while read id
do 
	if((i%$2==$3))
	then
  $id  
  fi # check the number1 number2
  i=$((i 1))
done
# for i in {0..9};do ( nohup bash ../submit.sh trim_galore.sh  10 $i 1>log.trim_galore.$i.txt 2>&1 &)  ;done

5 Hisat2比对

5.1 参考基因组下载

代码语言：bash复制
## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

# 进入到参考基因组目录
mkdir -p $HOME/database/GRCh38.111
cd $HOME/database/GRCh38.111

# 下载基因组序列axel  curl  
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz >dna.log &

# 下载转录组序列
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz >rna.log &

# 下载基因组注释文件
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gtf/homo_sapiens/Homo_sapiens.GRCh38.111.gtf.gz >gtf.log &

nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gff3/homo_sapiens/Homo_sapiens.GRCh38.111.gff3.gz >gff.log&

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
nohup gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz Homo_sapiens.GRCh38.cdna.all.fa.gz >unzip.log &

5.2 参考基因组建立索引

代码语言：bash复制
cd $HOME/database/GRCh38.111
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
vim Hisat2Index.sh

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

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

5.3 比对参考基因组

代码语言：bash复制
cd /home/data/t110558/project/Human-16-Asthma-Trans/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 -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

5.4 sam转bam

代码语言：bash复制
# 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
# 这个过程会输出大量中间文件
rm *.sam 
# 为bam文件建立索引
ls *.bam |xargs -@ 20 -i samtools index {}

同理，如果是样本数量太多了，就需要考虑并行并且控制提交样品数量

代码语言：bash复制
indexPrefix=/home/data/server/reference/index/hisat/mm10/genome 
ls -lh $indexPrefix*   
cat config|while read id;do
 if((i%$2==$3))
        then
  hisat2 -p 3  -x  $indexPrefix -1 ${id}*_R1_val_1.fq.gz   -2 ${id}*_combined_R2_val_2.fq.gz  -S ${id}.hisat.sam 
samtools sort -O bam -@ 3   -o ${id}.hisat.bam  ${id}.hisat.sam 
if [  ! -f  ${id}.hisat.bam    ]; then
		rm  ${id}.hisat.sam  
fi
samtools index ${id}.hisat.bam 
  fi # check the number1 number2
  i=$((i 1))
done 
  
# for i in {0..9};do ( nohup bash run_hisat2.sh t 10 $i 1>log.hisat2.$i.txt 2>&1 & )  ;done

6 FeatureCounts定量

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

代码语言：bash复制
cd $HOME/project/Human-16-Asthma-Trans/Expression/featureCounts

## 定义输入输出文件夹
gtf=/home/data/t110558/database/GRCh38.111/Homo_sapiens.GRCh38.111.gtf
inputdir=$HOME/project/Human-16-Asthma-Trans/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 &
代码语言：bash复制
# 对定量结果质控
multiqc all.id.txt.summary

# 得到表达矩阵
# 处理表头
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#/home/t_rna/project/Human-16-Asthma-Trans/Mapping/Hisat2//##g' |sed 's#.Hisat_aln.sorted.bam##g' >raw_counts.txt
#修改列名
sed -i '1s#/home/data/t110558/project/Human-16-Asthma-Trans/Mapping/hisat2//##g; 1s#.hisat.bam##g' raw_counts.txt
#查看表达矩阵
head -n 100 raw_counts.txt|column -t

得到的基因表达矩阵如图

数据挖掘

0 人点赞