从ncbi下载数据,GEO-GEO号/SRA-SRA号/SRA-SRP号,找到最终的SRP号,
在搜索界面send to file:Runinfo,内含https下载链接
或send to file:accession list,记下每一行值,补全链接为:
ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR(前三位)/SRR(全名)/SRR(全名).sra
下载目前使用idm和迅雷下载,prefetch下载太慢,ascp又不能用,后拷贝至分析端
0.环境准备(针对win10 ubuntu 18.04 LTS子系统)
linux清sudo密码(嫌烦)
sudo passwd -d username
apt换源(清华源)
cd /etc/apt/
sudo cp sources.list sources.list.bk
sudo vi sources.list
# dd删除一行,ggdG删除所有内容
清华源,有时会因为证书过期无法进行下去
# 默认注释了源码镜像以提高 apt update 速度,如有需要可自行取消注释
deb https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic main restricted universe multiverse
# deb-src https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic main restricted universe multiverse
deb https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-updates main restricted universe multiverse
# deb-src https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-updates main restricted universe multiverse
deb https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-backports main restricted universe multiverse
# deb-src https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-backports main restricted universe multiverse
deb https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-security main restricted universe multiverse
# deb-src https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-security main restricted universe multiverse
# 预发布软件源,不建议启用
# deb https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-proposed main restricted universe multiverse
# deb-src https://mirrors.tuna.tsinghua.edu.cn/ubuntu/ bionic-proposed main restricted universe multiverse
sudo apt-get update
sudo apt-get upgrade
安装图形界面
sudo apt install ubuntu-desktop
export DISPLAY=$(grep -m 1 nameserver /etc/resolv.conf | awk '{print $2}'):0.0 >> ~/.bashrc
export XDG_SESSION_TYPE=x11 >> ~/.bashrc
source ~/.bashrc
sudo service dbus restart
gnome-session
之后每次启动分为三步:
1.启动ubuntu20.04
2.启动Xlaunch,改动三处:选择one large window;display number设为0;勾上disable access control
3.gnome-session
安装中文语言包(随意)
安装miniconda3(4.6.14版本好使),手动激活环境变量(可选),更换清华镜像源
wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
chmod 777 Miniconda3-latest-Linux-x86_64.sh #给执行权限
bash Miniconda3-latest-Linux-x86_64.sh #运行
#echo 'export PATH="~/miniconda3/bin:$PATH"'>>~/.bashrc
#source .bashrc
#sudo chmod 777 -R miniconda3
#cd miniconda3/bin
#. ./activate
#一路默认
#初始化
conda config --set auto_activate_base false
conda config --set show_channel_urls yes
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
conda config --get channels
软件安装
conda install sra-tools htseq multiqc stringtie fastp macs2
sudo apt-get install fastqc bwa samtools bamtools trimmomatic tree pandoc unzip axel aria2 subread bowtie bowtie2 hisat2 (zsh,oh-my-zsh)
fastqc有时候运行有问题,最好去官网下载然后安装
trimmomatic如要使用官网下载即可
stringtie也去官网下载,git clone 然后 make
能用apt-get安装的最好不要使用conda安装,conda安装的有时运行起来有问题,比如说我的hisat2,bowtie2无法接samtools sort使用,很神奇,但用apt-get安装的就可以,可能和conda环境有问题,我还没找到解决办法
安装aspera
1.sra转为fastq
fasterq-dump -3 -e 8 -p -o 输出路径 sra文件
#(-e 指定线程数,-3<=>--split-3,-p 唠叨模式)
rename 's/.sra//' 输入文件(去除文件名中的.sra)
2.质检(现在我更喜欢fastp,速度快,长度可设置阈值,两端滑窗质量过滤,还有碱基校正)
fastqc -o 输出路径 -t 32 *.fq.gz
sudo su#运行multiqc需要提升权限
multiqc -t default -o 输出路径 -f 输入路径
java -jar trimmomatic-0.39.jar PE \
-threads 16 \
-trimlog trim.log \
input_forward.fq.gz input_reverse.fq.gz \
output_forward_paired.fq.gz output_forward_unpaired.fq.gz \
output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads \
SLIDINGWINDOW:5:20 \
LEADING:3 \
TRAILING:3 \
MINLEN:36
PE 表明处理的为双端数据,单端数据用SE
-thread 16 设置线程数为16
-phred33 设置碱基的质量格式(默认-phred64,自v0.32版本之后可自动识别是phred33还是phred64)
-trimlog trim.log 设置trimmommatic工具处理的日志文件为'trim.log',每两行为一对reads信息
input_forward_R1.fq.gz 输入的forward正向链R1对应的fastq文件
input_reverse_R2.fq.gz 输入的reverse反向链R2对应的fastq文件
output_forward_paired_R1.fq.gz 处理后输出的R1对应reads成对fastq文件
output_forward_unpaired_R1.fq.gz 处理后输出的R1对应reads不成对的fastq文件
output_reverse_paired_R2.fq.gz 处理后输出的R2对应reads成对fastq文件
output_reverse_unpaired_R2.fq.gz 处理后输出的R2对应reads不成对的fastq文件
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads 切除illumina接头参数设置。
SLIDINGWINDOW:5:20 设置滑动窗口阈值,以为5bp为窗口,这5bp碱基的平均质量值低于20,要进行切除
LEADING:3 设置从reads起始开始,去除质量低于阈值或为'N'的碱基,直到达到阈值不再去除,这里设置阈值为3
TRAILING:3 设置从reads末尾开始,去除质量低于阈值的碱基或为'N'的碱基直到达到阈值不再去除,这里设置阈值为3,这种方法是去除特定的illumina平台低质量区域(由于illumina会将低质量碱基标记为2),官方操作文档中建议使用 SLIDINGWINDOW 或 MAXINFO代替[这里未给出MAXINFO参数说明] )
MINLEN:36 设置read切除后的最短长度,这里设置长度至少为36bp,长度小于36bp的reads被去除
#fastp用法
fastp -w 32 \线程数
-5 \ 从5'端开始滑窗
-3 \ 从3'端开始滑窗
-w 4 \ 滑窗大小(长度)
-M 20 \ 滑窗平均最小质量值
-c \ 启用PE数据的碱基校正
-i fq1 \ fastq文件1
-o fq_1 \ clean后的fastq文件1
-I fq2 \ fastq文件2
-O fq_2 \clean后的fastq文件2
3.比对(参考基因组去hisat2官网下载,下载带“chr”的UCSC)
hiast2-extract-splice-sites.py Homo-sapiens.GRCh38.87.chr.gtf >hg38.ss
hiast2-extract-exons.py Homo-sapiens.GRCh38.87.chr.gtf >hg38.exon
hisat2-build -ss hg38.ss -exon hg38.exon 基因组文件位置 文件输出路径和名称
hisat2 --new-summary -x index位置 -1 1.fastq -2 2.fastq -S 文件输出路径和名称 2>x.log | samtools sort -@ 16 -o 1.bam(单线程比对效率高些)
4.排序并压缩为bam文件(可省略)
samtools view sam文件 -b > bam文件
samtools sort bam文件 -o 输出为排序后的bam文件(默认排序为染色体位置)
或执行下面:
samtools sort sam文件 -o 输出为排序后的bam文件(默认排序为染色体位置)
samtools index 输入为排序后的bam文件
5.计数(gencode的gtf和gff3都是带“chr”的)
htseq-count -f bam -r pos -s no -i gene_id 输入的bam文件 输入的gtf文件 > 输出的count文件
#(-f 输入的文件格式 -r 排序方式 -s 表示是否只比对单条链 -i GFF文件的一类属性,是计数单位,默认值是:gene_id)```
featureCounts -T 32 -F gtf -p -t exon -g gene_id -a 输入的gtf文件 -o 输出的count文件 输入的bam文件
#(-T 线程数 -F 注释文件类型 -p 指示结果为双端测序比对结果 -t 计数种类 -g 计数索引)