本文档记录GSE149638数据集中下载SRR11652578和SRR11652615原始数据
将其从SRR格式转换为fastq格式,通过fastqc质控,trim_galore过滤,Hisat2比对至基因组,FeatureCounts定量,最终得到基因表达的count矩阵的过程。
文件目录参考如下,每一个项目就在project中建立相应的文件夹
## 示例如下:
├── 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 # 差异分析
下载原始数据GSE149638
https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA629498&o=acc_s%3Aa
从96个数据集中下载最小的两个数据集作为测试
注:
“prefetch + fasterq-dump 的组合是从 SRA 访问中提取 FASTQ 文件的最快方法。预取工具将所有必要的文件下载到您的计算机上。如果下载不成功,可以多次调用 prefetch - 工具。它不会每次都从头开始;相反,它将从上次调用失败的位置开始。https://github.com/ncbi/sra-tools/wiki/08.-prefetch-and-fasterq-dump”
https://github.com/ncbi/sra-tools/wiki/02.-Installing-SRA-Toolkit
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
参考链接https://mp.weixin.qq.com/s/M8-ZuGgW2TQf2cKF4Mea3A
SRR_download_List.txt
## 记得激活环境
cat SRR_download_List.txt | while read id
do
echo prefetch ${id} -O ./
done > prefetch.command
less prefetch.command
nohup bash prefetch.command &
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 文件。以下是代码的逐行解释:
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
#样本数量少就可以直接提交到后台运行
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
#如果数量多,可以用这个脚本控制下并行数量
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 脚本内容如下所示;
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
## 参考基因组准备:注意参考基因组版本信息
# 下载,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 &
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 &
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
# 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 {}
同理,如果是样本数量太多了,就需要考虑并行并且控制提交样品数量
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
使用featureCounts对bam文件进行定量。
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 &
# 对定量结果质控
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
得到的基因表达矩阵如图
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。