whole-genome-sequencing Data Analysis 学习笔记4 变异位点注释数据库

linux下如何使用sftp命令
sftp> get /var/www/fuyatao/index.php /home/fuyatao/
这条语句将从远程主机的 /var/www/fuyatao/目录下将 index.php 下载到本地 /home/fuyatao/目录下。

sftp> put /home/fuyatao/downloads/Linuxgl.pdf /var/www/fuyatao/
这条语句将把本地 /home/fuyatao/downloads/目录下的 linuxgl.pdf文件上传至远程主机/var/www/fuyatao/ 目录下。

卸载软件：sudo apt-get autoremove –purge +软件名

变异数据库：对于找有意义的变异位点、疾病预测等方面有着重要作用

通常个体的全基因组测序数据可以挖掘到四百万个SNVs(跟参考基因组不一样的单碱基位点)
还有五十万的indels(insertions or deletions)
得到的数据通常是以vcf文件格式给出的

vcf格式
CHROM： 表示变异位点是在哪个contig 里call出来的，如果是人类全基因组的话那就是chr1…chr22，chrX,Y,M。

POS： 变异位点相对于参考基因组所在的位置，如果是indel，就是第一个碱基所在的位置。

ID： 如果call出来的SNP存在于dbSNP数据库里，就会显示相应的dbSNP里的rs编号。

REF和REF： 在这个变异位点处，参考基因组中所对应的碱基和研究对象基因组中所对应的碱基。

QUAL： 可以理解为所call出来的变异位点的质量值。Q=-10lgP，Q表示质量值；P表示这个位点发生错误的概率。因此，如果想把错误率从控制在90%以上，P的阈值就是1/10，那lg（1/10）=-1，Q=（-10）*（-1）=10。同理，当Q=20时，错误率就控制在了0.01。

FILTER： 理想情况下，QUAL这个值应该是用所有的错误模型算出来的，这个值就可以代表正确的变异位点了，但是事实是做不到的。因此，还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤，过滤完了之后，在FILTER一栏都会留下过滤记录，如果是通过了过滤标准，那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS，如果没有通过过滤，就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话，就说明没有进行过任何过滤。

#CHROM POS     ID        REF    ALT     QUAL FILTER INFO                              FORMAT      NA00001        NA00002        NA00003
20     14370   rs6054257 G      A       29   PASS   NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.

chr1 873762
GT: AD: DP: GQ: PL
0/1: 173,141: 282: 99: 255,0,255

其中最后面两列是相对应的，每一个tag对应一个或者一组值，如：

chr1：873762，
GT对应0/1；AD对应173,141；DP对应282；GQ对应99；PL对应255,0,255。

GT： 表示这个样本的基因型，对于一个二倍体生物，GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF一样；1表示表示跟ALT一样；2表示第二个ALT。当只有一个ALT 等位基因的时候，0/0表示纯和且跟REF一致；0/1表示杂合，两个allele一个是ALT一个是REF；1/1表示纯和且都为ALT；

0 - reference call
1 - alternative call 1
2 - alternative call 2

AD： 对应两个以逗号隔开的值，这两个值分别表示覆盖到REF和ALT碱基的reads数，相当于支持REF和支持ALT的测序深度。

DP： 覆盖到这个位点的总的reads数量，相当于这个位点的深度（并不是多有的reads数量，而是大概一定质量值要求的reads数）。

PL: 对应3个以逗号隔开的值，这三个值分别表示该位点基因型是0/0，0/1，1/1的没经过先验的标准化Phred-scaled似然值（L）。如果转换成支持该基因型概率（P）的话，由于L=-10lgP，那么P=10^（-L/10），因此，当L值为0时，P=10^0=1。因此，这个值越小，支持概率就越大，也就是说是这个基因型的可能性越大。

GQ： 表示最可能的基因型的质量值。表示的意义同QUAL。

举个例子说明一下：

chr1 899282 rs28548431 C T [CLIPPED]

GT: AD: DP: GQ: PL
0/1: 1,3: 4: 25.92: 103,0,26（该位点的基因型0/0，0/1，1/1）

在这个位点，GT=0/1，也就是说这个位点的基因型是C/T；
GQ=25.92，质量值并不算太高，可能是因为cover到这个位点的reads数太少，DP=4，也就是说只有4条reads支持这个地方的变异；
AD=1,3，也就是说支持REF的read有一条，支持ALT的有3条；
在PL里，这个位点基因型的不确定性就表现的更突出了，
0/1的PL值为0，虽然支持0/1的概率很高；
但是1/1的PL值只有26，也就是说还有10^(-2.6)=0.25%的可能性是1/1；但几乎不可能是0/0，因为支持0/0的概率只有10^(-10.3)=5*10-11

TCGA数据库下载脚本解析：
wget https://gdc-docs.nci.nih.gov/Data/Release_Notes/Manifests/GDC_open_MAFs_manifest.txt
for i in cut -f 2 GDC_open_MAFs_manifest.txt
do
echo $i
adress=echo $i |cut -d'.' -f 4

***cut -d “=” -f 2中的cut 是列提取命令，-d是判断文件是否存在的命令
-f表示取第一个字段的值。
如：echo “a/b/c” |cut -d ‘/’ -f 1，执行结果是a。执行过程：先按/分段，分段后结果是：第一个字段是a，第2个字段是b，第3个字段是c，-f就是取第几个字段。*

filename=echo $i |cut -f 2 |cut -d'.' -f 1-3,5-7
echo $adress$filename
wget -O “$filename” “”
done

for i in /*; do echo $i; find$i | wc -l; done;
这是一个for循环语句。
首先声明一个变量i ，将/*即根目录下的所有文件的名称全部赋值给刚刚声明的变量i。这一段在编程中叫做循环头。接下来。do 和done之间的部分是循环体。

echo命令的功能是在显示器上显示一段文字，一般起到一个提示的作用。
echo $i的作用是将根目录下所有文件显示出来。

wc -l 的作用是显示文件的行数。
find i是查找根