SMART-SEQ2特点:基因数量还行,细胞数量不太多。
这里是佳奥!让我们开始吧!
从下载文章数据开始。SRA是测序数据下载,Supplementary file是表达矩阵。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229
课程只对以下数据进行上游分析:
SRR6791441
SRR6791442
SRR6791443
SRR6791444
SRR6791445
SRR6791446
SRR6791447
SRR6791448
##从创建环境开始
conda create -n scrna
##添加到环境
export PATH="$PATH:/home/kaoku/biosoft/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin"
##建立srr.list并下载(下载目录ncbi/sra)
$ cat srr.list
SRR6791441
SRR6791442
SRR6791443
SRR6791444
SRR6791445
SRR6791446
SRR6791447
SRR6791448
cat srr.list | while read id; do ( prefetch $id & ); done
##数据下载完成
SRR6791441.sra SRR6791442.sra SRR6791443.sra SRR6791444.sra SRR6791445.sra SRR6791446.sra SRR6791447.sra SRR6791448.sra srr.list
##下载软件
conda install -y -c bioconda fastqc multiqc trim-galore subread hisat2
##SRR转fastq文件(单端测序),并新建raw_fq目录
ls /home/kaoku/project/scRNA/srr/*.sra |
while read id
do
fastq-dump -O ./ --gzip --split-3 $id &
done
##hisat2比对(mm10参考基因组)
index=/home/jianmingzeng/ reference/index/hisat/mm10/genome
ls raw_fq/*.gz |
while read id
do
hisat2 -p 10 -x $index -U $id -S ${id%%.*}.hisat.sam
done
##.sam转.bam
ls *.sam | while read id ; do (samtools sort -0 bam -@ 5 -o $(basename ${id} " .sam").bam ${id}); done
##构建index
ls *.bam | xargs -i samtools index {}
##count计数(需要.gtf文件)
gtf=/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf
featureCounts -T 5 -t exon -g gene_id -a $gtf -o all.id.txt *.bam 1>counts.id.1og 2>&1 &
##生成表达矩阵all.id.txt 与GEO数据库中的rawCounts结果相似
后续的分析我们就从解读作者提供的表达矩阵开始。
我们下一篇再见!