参考：转录组分析记录； 转录组入门和进阶
以下内容为转录组全部的上游分析，包括下载数据、比对、计数得到表达矩阵。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916
登录以上网址可直接进入目标GEO数据库（可直接更改最后的GSE号查找其他数据）
根据overall design，需要下载第9~15（rna数据）
image.png
image.png
image.png
image.png
image.png
进入SRA →all run,找到需要下载的数据，即SRR35899{56..62}
1.下载数据
(1)下载原始数据
①wget下载
#修改SRR号，逐一下载
wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra
#编写循环下载
vi wget.sh
for i in `seq 56 62`
do
nohup wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR35899${i}/SRR35899${i}.sra &
done
bash wget.sh
②prefetch 下载
a.下载安装Sratoolkit，并添加到环境变量
b.打开NCBI GEO ---> SRA ---> 任点一个进去 ---> All runs ---> Accession list --->得到需要的下载的所有SRA号
cat > srr.list
#打开Accession list，将里面的SRA号拷贝进srr.list
cat srr.list |while read id;do (nohup $prefetch $id &);done
(2)下载小鼠和人类参考基因组（fa格式）及注释文件（gtf格式）
以下演示的是小鼠基因组的下载，人类基因组类似
##注意都是从embl上下载
cd ../reference/genome/mm10/embl
for i in `seq 1 19` ;do wget ftp://ftp.ensembl.org/pub/release-96/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.${i}.fa.gz;done
wget ftp://ftp.ensembl.org/pub/release-96/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.X.fa.gz
wget ftp://ftp.ensembl.org/pub/release-96/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.Y.fa.gz
wget ftp://ftp.ensembl.org/pub/release-96/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.MT.fa.gz
将参考基因组合并为一个fa文件
cat Mus_musculus* > mm10.fa
#再下载gtf
从网站上直接下载：http://asia.ensembl.org/Mus_musculus/Info/Index
拷贝到../reference/gtf/mm10/embl
2.sra2fq
cd ../rna-seq/raw
ls ../rna-seq/sra/* |while read id;do (nohup fastq-dump --gzip --split-3 -O ./ $id & );done
3.校验
md5sum *fq.gz >md5sum
md5sum -c md5sum
4.对raw data质控（fastqc）
cd ../rna-seq/qc
##单个做质控
fastqc /data1/spider/liupiao/data/rna-seq/raw/SRR1039508.fq.gz -o ~/
##批量做质控
ls ~/raw/*gz |xargs fastqc -t 10
multifq ./
##将刚刚生成的qc文件综合在一起，跑个multiqc
##将生成的结果放在网页上查看
5.过滤（去掉接头和低质量序列）
cd ../rna-seq/raw
ls ../rna-seq/raw/*_1.fastq.gz >1
ls../rna-seq/raw/*_2.fastq.gz >2
paste 1 2 >configcat >qc.sh
source /data1/spider/miniconda3/bin/activate
dir='/data1/spider/liupiao/data/rna-seq/clean'
cat config |while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
##length可以根据qc看到reads的长度，做适当调整bash qc.sh config
6.比对（hisat2）
小鼠和人类的数据分开比对
①小鼠数据的比对（参考基因组为mm10）
##构建索引
cd /data1/spider/liupiao/data/index/hisat2/mm10
hisat2-build -p 4 /data1/spider/liupiao/data/reference/genome/mm10/embl/mm10.fa mm10
##开始比对
#注意小鼠数据是59~62。
for i in `seq 59 62` ; do
nohup hisat2 -x /data1/spider/liupiao/data/index/hisat2/mm10/embl/mm10 -1 ../SRR35899${i}_1.fastq.gz -2 ../SRR35899${i}_2.fastq.gz -S ../align/SRR35899${i}.sam &
done
②人类数据的比对（参考基因组为hg38）
7.sam2bam,对bam排序，生成索引
for i in `seq 59 62`
do
samtools view -S ../SRR35899${i}.sam -b > ./SRR35899${i}.bam
samtools sort ./SRR35899${i}.bam -o ./SRR35899${i}_sorted.bam
samtools index ./SRR35899${i}_sorted.bam
done
8.对bam进行质控（samtools flagstat ---> multiqc）
##构建索引ls *.bam |xargs -i samtools index {} ##生成.bai的索引文件ls *.bam |while read id ;do (samtools flagstat -@ 10 $id > $(basename ${id} ".bam").flagstat );done mkdir statmv *flagstat /data1/spider/liupiao/data/rna-seq/clean/test/statcd statmultiqc ./##在网站上看统计结果
9.计数
①htseq-count
文献中测序的数据是双末端PE-reads，htseq的计数需要进行按照reads名称进行排序
cd ../rna-seq/bam
conda /data1/spider/miniconda3/bin/activate
### samtools重新排序
for i in `seq 56 58`; do nohup samtools sort -@ 5 -n ../5_samtools/SRR35899${i}.bam -o SRR35899${i}_nsort.bam & done
###HTSEQ-count脚本
vim HTSEQ-count.sh
for i in `seq 56 58`
dohtseq-count -s no -r name -f bam SRR35899${i}_nsort.bam /data1/spider/liupiao/data/reference/gtf/hg38/embl/Homo_sapiens.GRCh38.95.chr.gtf > SRR35899${i}_matrix.count 2> SRR35899${i}.log
done
bash HTSEQ-count.sh
conda里面安装的HTSEQ-count总是出现库问题，无法运行
解决方案：pip
###退出conda环境
pip --version ##查看pip路径
###出现 pip 18.1 from /data1/spider/python3/lib/python2.7/site-packages/pip (python 2.7)
pip install HTSeq ##开始安装HTSeq及其依赖包
pip install HTSeq ##会提示已安装并显示安装的位置
###出现：Requirement already satisfied: HTSeq in /data1/spider/python3/lib/python2.7/site-packages (0.11.2)
cd /data1/spider/python3/lib/python2.7/site-packages
ls -lt
cd HTSeq
htseq-count ##出现用法等信息表明可以使用了
###因为是在Python下下载安装的，所有可以在任何路径中不输入绝对路径就能用了
##用前面写好的脚本
bash HTSEQ-count.sh
②featureCounts
##小鼠转录组
mm10_gtf=/data1/spider/liupiao/data/reference/gtf/mm10/embl/Mus_musculus.GRCm38.96.chr.gtf.gz
featureCounts -T 5 -p -t exon -g gene_id -a $mm10_gtf -o counts_mm10.txt *.bam
#人类转录组
hg38_gtf=/data1/spider/liupiao/data/reference/gtf/hg38/embl/Homo_sapiens.GRCh38.95.chr.gtf
featureCounts -T 5 -p -t exon -g gene_id -a $hg38_gtf -o counts_hg38.txt *.bam
网友评论