1.转录组 | 上游分析（数据下载，比对，计数）

2019-04-17 本文已影响23人 pomela

参考：转录组分析记录；转录组入门和进阶
以下内容为转录组全部的上游分析，包括下载数据、比对、计数得到表达矩阵。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916
登录以上网址可直接进入目标GEO数据库（可直接更改最后的GSE号查找其他数据）
根据overall design，需要下载第9~15（rna数据）

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