跳转至

2018

  • 分类于 肿瘤基因组
  • 需要 2 分钟阅读时间

初识Gene signature

什么是Gene signature

https://cancer.sanger.ac.uk/cosmic/signatures

Gene signature 的原理

人类的参考基因组序列是固定的，可以由此计算出任意三个连续的碱基序列的先的频率

  • 这样的三个连续的碱基ATG，TCG等共有96种
    • 密码子有64个，但是碱基组合有96（4 * 4 * 4 = 96）种
  • 参考基因组的这96中序列的出现的频率也是固定的，以柱状图显示出来，即使原始signature
    • 图为cosmic数据计算出的一种signature1

signature

  • 期吸烟的患者，体内基因发生变异，导致其测序统计后的signature与原始signature出现差异
  • 统计大量上述数据，得出的signature即可视为吸烟患者的signaure，cosmic数据给出了30种signature1-30
  • 测序获得一个新的signature-new,即可与signature1-30进行比较,从而解读该signature-new

Gene signature的用途

预后诊断

如何提取Gene signature

r语言，有两个两个包可以实现:

  • maftools
    • 这个包功能比较全面，不局限在signature上
  • deconstrucSigs
    • 仅仅实现signature作图，但是可以把signature的权重提取出来，示例如下

deconstrucSigs的安装

S
1
2
source("http://bioconductor.org/biocLite.R")
biocLite("")

输入数据格式

由病人的vcf文件中提取，整理为deconstrucSigs的格式

Text Only
1
2
3
4
5
6
7
Sample  chr      pos ref alt
1 chr1   905907   A   T
1 chr1  1192480   C   A
1 chr1  1854885   G   C
1 chr1  9713992   G   A
1 chr1 12908093   C   A
1 chr1 17257855   C   T
+ 注意decon.maf 中chr那一列为 chr1 ，不是 Chr1 + sep = "\t" + 多个病人的数据cat 在一起，用sampleID做区分

批量保存图片

其中 'XH01', 'XH11', 'XH12', 'XH19', 'XH28', 'XH31', 'XH33', 'XH34', 'XH36', 'XH38', 'XH39', 'XH40', 'XH43', 'XH44', 'XH45', 'XH47', 'XH49', 'XH51', 'XH52', 'XH54', 'XH56', 'XH58' 为我的数据中SampleID

S
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
suppressPackageStartupMessages(library("deconstructSigs"))
sigs.input <- mut.to.sigs.input(mut.ref="decon.maf",sample.id = "Sample", chr = "chr", pos = "pos", ref = "ref", alt = "alt")
l = list('XH01', 'XH11', 'XH12', 'XH19', 'XH28', 'XH31', 'XH33', 'XH34', 'XH36', 'XH38', 'XH39', 'XH40', 'XH43', 'XH44', 'XH45', 'XH47', 'XH49', 'XH51', 'XH52', 'XH54', 'XH56', 'XH58')
pdf("plot.cosmic.pdf")
for(i in l){
test = whichSignatures(tumor.ref = sigs.input,signatures.ref = signatures.cosmic,sample.id=i,signature.cutoff=0,contexts.needed = TRUE,tri.counts.method = 'default')
plotSignatures(test)
makePie(test)
}
dev.off()
+ 该方法有个缺陷，图片输出的title字体太大，当signature较多时显示不完全

批量输出各signature 的weights
S
1
2
3
4
5
for(i in l){
    test <- whichSignatures(tumor.ref = sigs.input,signatures.ref = signatures.cosmic,sample.id=i,signature.cutoff=0,contexts.needed = TRUE,tri.counts.method = 'default')
    res <- test$weights
    write.table(res,"cosmic.weight.txt",sep="\t",col.names=F,row.names=T,quote=F,append=T)
}
  • 可以将 signatures.ref = signatures.cosmic 改为 signatures.ref = signatures.signatures.nature2013 ,这是两个不同来源的数据集算出的signature，我们用计算出的signature 与该数据的signature进行比较，为我们的数据做注释
    • 得出的结果是：signature_sample1 = signature1 * weight1 + signature2 * weight2 + ··· + signature30 * weight30
    • 通过权重的比值对我们算出的signature进行注释
  • 分类于 生物信息学
  • 需要 5 分钟阅读时间

生物信息学常用的Linux命令

1、 前言

东西有点杂乱，还需要慢慢整理。其中的 sed awk perl 等命令我觉得都可以单独分出来讲了

2、 构建blast数据库

Bash
1
2
3
4
makeblastdb -in ../ATaa.fa -dbtype prot -title ATaa-parse_seqids -out ATdb   #可用的蛋白库
makeblastdb -in ../ATaa.fa -dbtype nucl -title ATaa-parse_seqids -out ATdb  #可用的核酸白库
blastp -db /home/hanyapeng/prot/DB/Mtrdb -query At.fa -outfmt 6 -evalue 1e-20 -out Mtr.blast -num_threads 8
blastn -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8

3、 碱基/蛋白序列比对

Text Only
1
2
3
muscle -in AtLj.fa -out AtLj.phy  # -out，输出默认的fasta格式
muscle -in At.fasta -phyiout At.phy  # -phyiout，则输出sequential phylip格式
mafft aa.fa > aa.fas

4、 画树

Bash
1
2
3
nohup raxmlHPC-PTHREADS -s ALO.phy -n ALO -m PROTCATJTT -f a -N 100 -p 12345 -x 1234567 -T 8 &  # RaxML树
nohup ../software/PhyML_3.0_linux64 -i AGO.phy -d aa -b 0  -f e -s BEST -c 4 -m VT -v e -a e -o lr --no_memory_check &  # Phyml 树
../software/FastTree TCP123.phy  > TCP.tree   #fasttree 画树

5、 提取CDS/pep序列

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
nohup hmmsearch -E 1e-5 TCP.hmm  /home/hanyapeng/prot/Gma.prot > Gma.hmm &
less -S Gma.hmm |awk '{if ($1~/e-/) print $9}' > list     #if第一行匹配到e-，则调出第9行#  

nohup hmmsearch --domtblout hmmresult.txt -E 1e-5 TCP.hmm  /home/hanyapeng/prot/Gma.prot  &
less -S hmmresult.txt |awk '{ if ($1~/e-/) print$1"\t"$20"\t"$21}'|less -S > GmaTCPdomain.list    # 提取第1.20.21列为坐标

perl /home/lixiangyong/software/getseq.pl Osa.list /home/hanyapeng/prot/Osa.prot > Osa.fa
perl /home/lixiangyong/software/get_domain.pl  /home/hanyapeng/prot/Osa.prot  OsaTCPdomain.list > OsaTCPdomain.fa  #提取domain

perl /home/lixiangyong/software/getseq.pl Chr13 Gmax275.fa > Chr13.fa  
perl /home/hanyapeng/comm/getpromoter.pl /data/plant_genome/Glyma2.0/Gmax275.fa Q1.txt > 1  # 提取启动子序列 Q1为启动子坐标
cat AT.fa Lj.fa Osa.fa > ATLjOsa.fa

6、 修改环境变量

Text Only
1
2
3
4
5
6
vi .bash_profile     # 或则 .bash_login   .profile   .bashrc 
vi .bash_profile     # 添加环境变量 PATH=$PATH:$HOME/bin
source .bash_profile  #   立刻执行更改

nano .bashrc             #和vi类似的编辑器
source .bashrc

7、 prottest的使用

Text Only
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
nohup java -jar prottest-3.4.jar -i AGO.phy -all-matrices -all-distributions  > a.txt  &  # 此.phy 文件为比对过的氨基酸序列文件

nohup java -jar prottest-3.4.jar -i namefail.phy -o namefail.txt -S 1 -AIC -I -G -IG -F -ncat 8 -all &

nohup java -jar prottest-3.4.jar -i AGOdomain.phy -o AGOdomain.txt -all-matrices -all-distributions &

nohup java -jar jModelTest.jar -d example-data/aP6.fas -f -i -g 4 -s 11 -AIC -a >a.txt  &  此.fas文件为核苷酸序列文件

PhyML-3.1_linux64 -i aP6.phy -d nt -n 1 -b 0 --run_id F81+I -m 000000 -f m -v e -c 1 --no_memory_check -o tlr -s BEST

PhyML-3.1_linux64 -i aP6.phy -q -d aa -m JTT -c 4 -a e

8、 shell 脚本知识

Bash
1
!= 不等于,如:if [ "$a" != "$b" ]

9、 awk的使用

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
less Y2.blast |awk ' {if ($11 ~ /0\.0/) print $0}' > 1      # 每一行中，如果第11列=0.0，则输出全部  “~”是 =  ，$11 调取第11列 ， \ 反义符 
awk '{line[NR]=$0}END{for(i=NR;i>0;i--)print line[i]}' sed.txt  # 把sed.txt 中的内容按行倒着输出
less -S Gma3.fa |awk '{if ($1~/>/) print $head}' > a
awk '!seen[$0]++' <filename>    # awk 去除重复的行
less Step1.2.sh | grep "=" | cut -d "=" -f 1 | awk '{printf $quot" ";$0}' #将一列输出为一行
awk '{ print $(NF-1)}' file     # 出输出倒数第二列

less file |  awk  '{split($1,a,":");print a[1]"-"a[3]"\t"$5}' > our.w.d    # 对awk截取字符串，file中的数据格式为:`ACAN:ENST00000561243:p.D1390E The mutant peptide: APGVEEISGLP is found in wild type sequence`， [参考](https://www.cnblogs.com/sunada2005/p/3493941.html)

awk 'BEGIN{total=0}{total+=$1}END{print total}' file   # 对文件第一列求和

10、 sed的使用

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
sed -i 's/Gma/Glyma/g' Glyma.fa
sed -i '/^$/d' 2_domain.txt  # 删除文件空行
sed -i '/#/d' AtTCP.txt      # 删除有\#的行
sed -i s/^M// AAt.list.txt   # 此处^M必需手动输入！！
sed -i /^$/d file     #删除空行   d 指删除含有空行符的这一行
sed -i s/tr.*//g LjTCP.fa  #  去除无用的注释    :%s/|.*.\[//g  #去除从|开始，到[结束内的所有字符，包括[#

sed 's/^/aaa/g' filename.txt > outputfile.txt  # 在文件中的每一行之前加上一个字符串，比如：aaa    行尾将^是 $ 
sed -i/\>Chr17/d Chr17.fa   # 删除Chr17.fa 中的Chr字符
sed -i "s/zhangsan/lisi/g" `grep zhangsan -rl modules` #批量替换，讲文件夹modules中所有文件执行sed -i       "s/zhangsan/lisi/g" 该符号`不是点，是Esc下面的符号 
sed -i '1iNew Top Line' aa.txt    #  sed 在句首添加一行 New Top Line

11、 linux 命令

Bash
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
date   #  查看系统日期
cal 1 2016   # 查看2016年1月的日历
df     # 查看磁盘驱动器的可用空间  
free   # 显示可用内存
.     # 工作目录     
..   # 工作目录的父目录   如：  cd ..   返回上一级    cd  ./software/  
file Gma1.fa   # 显示  Gma1.fa: ASCII text  该文件为文本文档ASCII
alias foo='cd /usr; ls; cd -'   # 自己创建命令   unalias foo 删除
gedit Gma1.fa   # 图形编辑器编辑Gma1.fa 

rm -- -foo   # 要删除第一个字符为‘-’的文件 (例如‘-foo’)
rm ./-foo   # 要删除第一个字符为‘-’的文件 (例如‘-foo’)

chmod 755 lixiangyong/  # 更改目录权限为可读

wc -l filename # 就是查看文件里有多少行
wc -w filename # 看文件里有多少个word。
wc -L filename # 文件里最长的那一行是多少个字
wc -c  filename   # 统计字节数

cat Y.list | sort -u > YRedundancy.list    # 去除文件Y.list中的重复的行命名为YR.list
less -S AGO.fa |grep ">" |wc -l  # 检索>有多少个  同  grep -c ">" AGO.fa
for tar in *.tar.gz;  do tar xvf $tar; done    #批量解压文件#
split -l 10 splitfile.txt D  # 将文件按每10行分割，新文件命名为 Daa Dab#
rename linux unix linux*    # 文件批量  
find . -name '*.txt' | xargs perl -pi -e 's/aaa/bbb/g'  # 遍历整个文件夹，把所有txt文件中的aaa替换成bbb
mkdir $(date +%F)  #创建一个以日期为名的文件夹
touch {1..9}.txt   # 批量创建文件

ls -F | grep '/$'         #  查看当前目录下的文件夹 
for i in `ls -F | grep '/$'`;do sudo du -sh $i;done   #查看当前目录下，每个子目录及其子目录的总大小

rename "s/AA/aa/" *
rename .fas .fa *.fas    # CentOS 下可用

iostat -d -k 1 10可以看到每个盘的读写速度

find /home/lxy -type f -size +100M    #列出文件夹 /home/lxy 下大于100M的文件
find . -type l -exec rm  {} \;  #列出当前路径 .  下的link类型的文件(type
 l) 并删除

12、 单行 perl

Text Only
1
perl -e "print lc($a)"    # Perl 大小写转换

13、 R 命令行

Text Only
1
Rscript /home/muscle/test.r  # linux 上运行r脚本

14、 kaks计算

Text Only
1
2
perl pal2nal.perl aligned.protein.fas original.cds.fas -output fasta > aligned.cds.fas
perl ../software/pal2nal.pl 2.fas 2.cds -output paml > 2.nuc    #kaks计算

15、 gff和gtf之间的格式转化

Text Only
1
2
gffread my.gff3 -T -o my.gtf # gff2gtf
gffread merged.gtf -o- > merged.gff3  # gtf2gff

16、 vim 的使用

Text Only
1
2
3
4
5
6
7
8
9
:12,23s/aa/bb/g     # 将从12 行到23 行中出现的所有包含aa 的字符串中的aa 替换为bb
:s/\<aa\>/bb/g      # 将光标所在行出现的所有aa 替换为bb, 仅替换aa 这个单词
:12,23s/^/#/        # 将从12 行到23 行的行首加入# 字符

:10,15s#^#\#\##g          #  vim的 V模式下，下注释x  10 -15  行，##注释

Esc shift+ : wq   enter   =  Esc shift+ : x  enter       #保存并退出
Esc shift+ : q！  enter                                  #不保存强制退出 
vi foo.txt ls-output.txt    # 同时vi编辑多个文件  ：n   ：N   ：buffer 2 来回切换     也可以先打开一个，在用：e ls-output.txt

17、 数据存放

/data/public/hanyapeng/Gm/Gma.collinear.groups # 大豆中所有旁系同源基因
/data/public/wanglei/inparanoid4.0/Soybeanparalog2.txt # 大豆中所有旁系同源基因

推荐连接

  • 分类于 博客
  • 需要 1 分钟阅读时间

搭建记录

  1. 在github账号下创建 biolxy.github.io
  2. git clone 该项目到本地
  3. git clone 作者博客到本地，填补博客必要的样式文件,删除_post文件夹下文件，保留一个做模板，命名最好为全英文字符
    Bash
    1
2
    git clone https://github.com/Gaohaoyang/gaohaoyang.github.io.git 
mv gaohaoyang.github.io/* biolxy.github.io
  4. 修改如下等文件替换为自己的信息
    CNAME _config.yml favicon.ico index.html page/3collections.md age/4about.md
  5. 该主题更详细的信息需要去该主题作者的github查看地址
  6. push到github，打开 https://biolxy.github.io/ 即可看到效果
  7. 当然你也可以
  8. 先参考主题作者的博客地址，先 安装 ruby 和 jekyll 环境 ，博客里讲的都有，我就不赘述了
  9. 切换到biolxy.github.io 文件夹下，注意不是_post 文件夹
  10. 执行 jekyll s
  11. 此时打开浏览器，输入http://127.0.0.1:4000/ ，你就可以在本地预览到上传github后的效果了
  12. 为了上传图片， 你需要注册七牛云，建立空间，做图床（我写的简单，你操作的时候不懂的就百度）
  13. 七牛云账户可以绑定MPic软件，绑定后可以方便的上传图片及插入图片链接。

更新

  • 七牛云不能用了，原来部分图床可能挂掉了 解决办法：
  • blog可以先发布在CSDN上，github上的图片直接连接CSDN的博客中的图片地址（为了扩大影响，都是一文多平台发布）
  • 分类于 生活札记
  • 需要 1 分钟阅读时间

一个test博客

1、 段落标题

一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧，一个段落，好像有点短啊，那就cp几遍吧
嗯，够了吧
得瑟

2、 列表测试

  • 冒烟：一个名不见经传的人。
  • 可爱：一个天生丽质的人。

3、 表格测试

书籍 技能
Python高性能编程 python

4、 一段代码

Python
1
2
3
4
5
6
7
@requires_authorization
class SomeClass:
    pass

if __name__ == '__main__':
    # A comment
    print 'hello world'

5、 公式

\[a + b = c$$ $\overbrace{a+b+\dots+n}^{m个}$ 行中 $$ f(x) = ax + b \]