原始数据来源于这篇文章https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50177
这篇文章的数据适中，不仅可以用来做RNA-seq，后面我们还可以用来做CHIP-seq
1.用Aspera下载原始数据并存放到Sra文件夹里面
~/.aspera/connect/bin/ascp -T -i /home/qiujunhui/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m
anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957677/SRR957677.sra ./
~/.aspera/connect/bin/ascp -T -i /home/qiujunhui/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m
anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957678/SRR957678.sra ./
~/.aspera/connect/bin/ascp -T -i /home/qiujunhui/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m
anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957679/SRR957679.sra ./
~/.aspera/connect/bin/ascp -T -i /home/qiujunhui/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m
anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR957/SRR957680/SRR957680.sra ./
如果下载出现防火墙问题，则需要手动打开端口33001
iptables -I INPUT -p udp --dport 33001 -j ACCEPT
iptables -I OUTPUT -p udp --dport 33001 -j ACCEPT
2.下载hisat2索引文件
用hisat2进行比对，我们需要建立索引文件，这里我们的索引文件是人类基因组hg19的索引文件，可以去官网下载
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
mkdir hg19
mv hg19.tar.gz hg19
tar -zvxf hg19.tar.gz #将索引文件解压存放在hg19文件夹里面
3.下载人类基因组注释文件
在用htseq-count计算基因表达量的时候，我们需要用到参考基因组的注释文件,可以在gencode数据库下载
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/GRCh37_mapping/gencode.v29lift37.annotation.gtf.gz
mv GRCh37_mapping/gencode.v29lift37.annotation.gtf.gz hg19
tar -zvxf GRCh37_mapping/gencode.v29lift37.annotation.gtf.gz #将基因组注释文件保存至hg19文件夹并解压
4.用脚本实现批量解压>比对>计算基因表达量一体化
为了简化操作我自己写了个脚本
#!/bin/bash
#把sra文件都存放在Sra文件夹里面
#参考基因组和注释文件都在hg19文件夹里面
#qiujunhui 1801963472@qq.com
mkdir Fastq_out #创建一个存放SRA-toolkit解压sra文件的输出文件夹
mkdir Sam_out #创建一个存放用hisat2比对的输出文件夹
mkdir Sorted.bam #创建一个存放用samtools将sam文件排序的生成的bam文件的文件夹
mkdir Counts #创建一个存放用HTseq-count计算基因表达量的输出文件夹
#批量解压文件
for i in ~/Sra/*sra
do
echo $i
#判断sra文件时单端测序还是双末端测序
num=$(fastq-dump -X 1 --split-spot -Z $i | wc -l | grep [0-9])
if [ $num -eq 4 ];then
echo "$i是单端测序"
fastq-dump $i
mv ~/Sra/*fastq ~/Fastq_out
#用hisat2进行比对
for x in ~/Fastq_out/*fastq
do
echo $x
a=$(echo $x | cut -d "." -f1)
hisat2 -p 5 -x ~/hg19/genome -U $x -S $a.sam
mv ~/Fastq_out/*sam ~/Sam_out
done
else
echo "$i是双端测序！"
fastq-dump --split-files $i
mkdir 1Fastq_out
mkdir 2Fastq_out
mv ~/Sra/*_1.fastq ~/1Fastq_out
mv ~/Sra/*_2.fastq ~/2Fastq.out
#用hisat2进行比对
for j in ~/1Fastq_out/*_1.fastq
do
for h in ~/2Fastq_out/*_2.fastq
do
b=$(echo $j | cut -d "_" -f1)
c=$(echo $h | cut -d "_" -f1)
if [ $b = $c ];then
hisat2 -p 5 -x ~/hg19/genome -U -1 $j -2 $h -S $b.sam
mv ~/2Fastq_out/*.sam ~/Sam_out
fi
done
done
fi
done
#用samtools进行排序
for y in ~/Sam_out/*sam
do
echo $y
d=$(echo $y | cut -d "." -f1)
samtools sort -n -@ 5 -o $d.bam $y
mv ~/Sam_out/*bam ~/Sorted.bam
done
#用HTseq-count进行基因表达量计算
for z in ~/Sorted.bam/*bam
do
echo $z
e=$(echo $z | cut -d "." -f1)
htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty $z ~/hg19/gencode.v29lift37.annotation.gtf>$ecounts.txt
mv ~/Sorted_bam/*.txt ~/Counts
done
#把不要的文件夹删除只留下counts文件夹
rm -rf ~/Sra
rm -rf ~/Fastq_out
rm -rf ~/1Fastq_out
rm -rf ~/2Fastq_out
rm -rf ~/Sam_out
rm -rf ~/Sorted.bam
计算出来的基因表达量存放在Counts文件夹里面
我们利用虚拟机和电脑的交互软件将文件保存到桌面，便于下一步用DESeq2进行基因表达差异分析
网友评论