序列的下载与指控
需要的序列
这一步主要是从各大数据库下载别人做好的高通量测序，可能一个样本会测好多次，就要下好多个样本
这次我下的是董老师要求的GSE99531（具体这个基因的研究回头再看）链接如下GSE99531
他的高通量测序是SRP108393，可以看到一共是测了37个样本
然后Send results to Run selector一下，在里面有个Accession List可以把所有序列的名字导出来。
后面是下载了
今天是20200525，没错这个下载序列一下子折腾我三天，前天晚上开始下了，期间试了不少乱七八糟的方法。
一共是用了几个工具，下面一一说一下吧。
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文件
这儿反正常规手段不行了x、
附：找人帮忙
还问了一个加拿带的朋友，想让他帮我下的，问了一下速度，随随便便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后再批量质控。
