序列的下载与指控

需要的序列

这一步主要是从各大数据库下载别人做好的高通量测序，可能一个样本会测好多次，就要下好多个样本

这次我下的是董老师要求的GSE99531（具体这个基因的研究回头再看）链接如下GSE99531

他的高通量测序是SRP108393，可以看到一共是测了37个样本

然后Send results to Run selector一下，在里面有个Accession List可以把所有序列的名字导出来。

后面是下载了

今天是20200525，没错这个下载序列一下子折腾我三天，前天晚上开始下了，期间试了不少乱七八糟的方法。折腾到今天，不得不说心阶的VPN质量真的好。

一共是用了几个工具，下面一一说一下吧。

sra-tool

这个应该是NCBI官方的一个东西吧。可以去官网看一下怎么用。SRAtool网页参考

这儿用的主要命令就是一个preftech

如果是单个测序文件的话，就是prefetch SRR5635554，多个测序文件可能有点麻烦，首先搞到上一个部分弄到的序列名字导出的SRR_Acc_List.txt，然后cd到这个目录下，之后写一个小循环，这里我抄的生信技能树视频里的。

cat SRR_Acc_List.txt |while read id ;do prefetch $id;done

这句就是，cat读取这个txt里的数据，然后用read逐行读取，并将读到的值放到id中，然后prefetch一下每个id。

这儿写了一个循环，语法应该是while [command] ; do [command] ; done

这个是我最后使用的一个命令，应该是最原始的下载NCB数据的操作，但是因为服务器在米国，所以速度非常感人。因为用的虚拟机，后来又折腾了一下翻墙才把速度提上来。

Aspera

• 一开始prefetch速度非常慢，然后我就考虑有没有什么快的方法，然后我去生信技能树求助了一下，虽然最后没找到结果，但是看到了这个工具。所以抱着试一试的心态去看了看

这个可以拿去做参考，基本上按照这个流程就可以：安装Aspera Connect工具下载sra数据

一般来说，NCBI的sra文件在ftp网站里，前面的路径都是一样的/sra/sra-instant/reads/ByRun/sra/SRR/SRR数据前3位/SRR全长/SRR全名.sra

所以我这个路径就是/sra/sra-instant/reads/ByRun/sra/SRR/SRR563/SRR5635554/SRR5635554.sra

然后这个流程有个问题，就是33001端口，基本上路由器都没开，我估计意思就是像ftp那样要把这个端口开起来。如果是在实体机上应该做个转发就行，但是虚拟机就，根本不会搞，尝试着把下面两行打上去但是没用

sudo iptables -I INPUT -p tcp --dport 33001 -j ACCEPT
sudo iptables -I OUTPUT -p tcp --dport 33001 -j ACCEPT

实体机也许可以用？

• 还有一个是直接上NCBI的ftp，但是好像没权限，要有aspera权限才可以上，但是这个插件对于Win和Mac只支持3.6.6以前的版本。我找了一圈没找到这么早之前的插件所以放弃了。

wget

网上说好像是可以的，但是我基本上等于没成功（可能太久之前的情报了），就放弃了，贴个链接给有缘人尝试吧。使用wget下载NCBI GEO数据库的sra文件

给实体机翻墙给虚拟机上

这个是prefetch最后的解决办法，找了个节点把问题基本解决了。vm虚拟机下如何实现ubuntu代理上网

ubuntu界面如何翻墙

因为之前一直不是很懂代理的机制，所以只是会用，这次换了ubuntu20之后根本用不起来了，所以后面换回了ubuntu18来试试。换v2ray之后感觉良好，最后用的Qv2ray。https://qv2ray.github.io/ 注意这个软件要下v2ray的框架的。

附：找人帮忙

还问了一个加拿带的朋友，想让他帮我下的，问了一下速度，随随便便4~5M/S，呜呜呜呜呜羡慕死了

小插曲

中间下到一半磁盘容量不够了呜呜呜，后来没办法只能关了扩了一次磁盘，真没想到这个数据能有50G。（好不容易遇到一次高速，舍不得断开）

sra转fastq

nohup ls *.sra|while read id; do (fastq-dump --gzip --split-3 -O /home/data/vip14/dick/a549/2_fq/ $id &); done > ../2_fq/nohup.log 2>&1 & 

注：这里必须要有个括号，不然“&”和分号在一起会报错。（原因不明）
解析一下：我所有sra都在circ文件夹下，就搞了个循环对所有sra都做fastq-dump这个命令，--gzip把文件压缩，--split-e出两条fastq，-O导出到文件夹。因为是双向测序，所以会出两个文件。

zless可以将内容显示。

质控

基础语句

fastqc SRR5635537_1.fastq.gz -O ./

没有-O可能会导出到一个没权限的文件夹？

多个一起处理就是

ls *gz |xargs fastqc -t 20 -O ./

-t是可以指示多线程。我电脑好像是32线程的？我跑了之后好像是跑满4个线程。（我四个文件）

如果质控的fastq太多了，生成太多html，会导致没法一个个看，所以要multiqc，语句把所有质控整合到一个表上

前提：所有fastqc都已经跑完了

multiqc ./

去adapter

因为测序序列会有很多adapter，为了去除adapter，要用一个软件来完成这一步。

去完adapter后再质控。

使用代码如下，各个参数懒得翻译了。

trim_galore -q 25  --phred33 --length 25 -e 0.1 --stringency 3 --paired -o dir fq1 fq2

-q/--quality <INT>：Trim low-quality ends from reads in addition to adapter removal.

--phred33：Instructs Cutadapt to use ASCII+33 quality scores as Phred scores(Sanger/Illumina 1.9+ encoding) for quality trimming. Default: ON.

--length <INT>：Discard reads that became shorter than length INT because of either quality or adapter trimming. A value of '0' effectively disables this behaviour. Default: 20 bp.For paired-end files, both reads of a read-pair need to be longer than <INT> bp to be printed out to validated paired-end files (see option --paired). If only one read became too short there is the possibility of keeping such unpaired single-end reads (see --retain_unpaired). Default pair-cutoff: 20 bp.

-e <ERROR RATE>：Maximum allowed error rate (no. of errors divided by the length of the matching region) (default: 0.1)

批量：因为我数据量小所以我，没用过这串？

ls ./raw/*_1.fastq.gz >1
ls ./raw/*_2.fastq.gz >2
paste 1 2 > config
dir='./clean'
cat $config_file |while read id
do
    arr=($id)
    fq1=${arr[0]}
    fq2=${arr[1]}
trim_galore -q 25  --phred33 --length 25 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2
done

跑完trim往回找一步，重新做质控

注：我觉得可以拿一两个出来坐下质控看下有没有adapter，然后一起trim后再批量质控。