加载中…
个人资料
邵先生
邵先生
  • 博客等级:
  • 博客积分:0
  • 博客访问:89,989
  • 关注人气:35
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
分类
搜博主文章
访客
加载中…
评论
加载中…
博文
(2020-11-25 15:31)
你要的群体样本PCA作图我已经帮你写好了,请直接用!!
嗯,就是方便,快捷!!

rawdata=read.table('final_matrix.pca',header=T)
##读取文件,行为基因,列为样本
row.names(rawdata) <- rawdata $gene
tmp <- t(rawdata[,c(-1)])
#第一列的基因名去掉
cleandata <- tmp[,colSums(tmp != 0) > 0]
#去掉行全为o的行
data.pca <-prcomp(cleandata, center= T,scale.= F)
#用prcomp做PCA,输出特征向量,

png('pca-0.png')
screeplot(data.pca,type='line',main='碎石图',lwd=2)
dev.off()
#画个碎石图看看

png('pca-1.png')
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2020-09-03 08:09)

说明:本程序主要是在小窗口下能高清显示图片能特意定制的参数,里面涉及的数值可以根据具体项目设置


阅读  ┆ 评论  ┆ 转载 ┆ 收藏 

第一步:格式转换

vcftols --vcf sample.vcf --plink --out sample

plink  --file sample --chr-set 20 --noweb --make-bed --out  sample

第二步:Q矩阵和画图

for i in 2 3 4 5 6 7 8 ;do python  structure.py -K $i --input sample --output  sample  --full --seed=100;done

1For i in 2 3 4 5 6 7 8;python distruct.py -K $i --input=sample  --title='' --output=$i.structure.svg;done

2http://omicsspeaks.com/strplot2/网页版,但是输入必须是csv格式

 

阅读  ┆ 评论  ┆ 转载 ┆ 收藏 

Gene-Based Association (GBA) analysis是基于全基因组关联分析的一种目前GWAS可以分为基于突变(SNP,CNV,SV),基因和通路三种。其中基于突变的主要是单核苷酸多态(SNP),目前,可以使用的软件比较多,比如TASSEL,GAPIT,PLINK等。从近几年发的文章来看,和人相关的GWAS分析一般都是使用PLINK,而动植物的除了TASSEL,GAPIT外,还有EMMAX,GEMMA,FastLMM等。支持gene-based关联分析研究的工具相对来说偏少,本人翻了google后发现有两个软件比较适用,分别是GCTA和VEGAS2,具体使用方法如下:

1、GCTA

参考资料:https://gcta.freeforums.net/board/2/gcta-user-manual

下载地址:https://cnsgenomics.com/software/gcta/#Download

第一步:数据格式转换

plink --allow-extra-chr --vcf test.vcf --make-bed --out test

plink --file test --make-bed  --out test --recode --allow-extra-chr

第二步:fastBAT

gcta64 --bfile test --maf 0.05 --fastBAT assoc.txt --fastBAT-gene-list gene_list.txt --out test --thread-num 10

其中:

(1)assoc.txt:

1:

阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2019-11-11 08:42)
文章摘自:https://zhuanlan.zhihu.com/p/87408242
https://github.com/stat-lab/EvalSVcallers


Option 3: Install using a Python Virtual Environment
    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
    (2019-07-02 15:02)

    1、背景 

    超敏位点和基因表达有关,并且超敏位点反应了染色质的可及性。也就可以反推出“可及性”的染色质结构区域可能与基因表达调控相关。使用了超敏Tn5转座酶切割染色质的开放区域,并且加上接头(adapter)进行高通量测序。一共需要三种酶,能切割出单个核小体的MNase, 能识别超敏位点的DNase ATAC-Seq所需要的

    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
    (2019-03-28 15:41)

    值得注意的是: All reads with a mapping quality < 70 were removed prior to calling.

    其它策略包括:

    • Based on our analysis of replicates, SNVs with MuTect quality scores <6.95 were removed.
    • We removed those variants that overlapped with repetitive regions
    • Fisher’s exact test was used to identify variants exhibiting read direction bias
    • SNVs present at VAFs smaller than 0.1 or at loci covered by fewer than 10 reads were removed, unless they were also present and confirmed somatic in the Catalogue of Somatic Mutations in Cancer (COSMIC).
    • 删除那些在千人基因组计划的任意人群(AMR, ASN, AFR) 里面频率大于1%的变异位点。
    • We used the normal samples in our data set (normal pool) to control for both sequencing noise and germline variants, and removed any SNV observed in the normal pool (at a VAF of at least 0.1).

    主要是区分recurrent和inact

    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
    邵先生有话说:在工作中难免领导会给你一堆数据,然后让你去分析,很多情况下不清楚是用什么仪器测序的,质量体系是多少,这时候就是抓心挠肺的时候,不用急,以下脚本帮你解决这个问题(这个脚本我也不清楚哪来的,如有雷同可以留言申明)

    less $1 | head -n 1000 | awk '{if(NR%4==0) printf('%s',$0);}'| od -A n -t u1 -v | awk 'BEGIN{min=100;max=0;} \
    {for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i
    {if(max<=126 && min<59) print 'Phred33'; \
    else if(max>73 && min>=64) print 'Phred64'; \
    else if(min>=59 && min<64 && max>73) print 'Solexa64'; \
    else print 'Unknown score encoding'; \
    print '( ' min ', ' max, ')';}'

    eg:
    输入:
    sh  fq_qual_type.sh
    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
    (2019-03-18 16:56)
    邵先生有话说xargs的功能可以用其它组合命令代替,但是如果使用xargs会使命令行更简洁,是计算更简单,如果你还没有看到这个命令的话,就一起来学习吧。

    NAME
           xargs - build and execute command lines from standard input

    SYNOPSIS
           xargs  [-0prtx]  [-E  eof-str] [-e[eof-str]] [--eof[=eof-str]] [--null]
           [-d delimiter] [--delimiter delimiter]  [-I  replace-str]  [-i[replace-
           str]]    [--replace[=replace-str]]   [-l[max-lines]]   [-L   max-lines]
           [--max-lines[=max-lines]] [-n max-args] [--max-args=max-args] [-s  max-
           chars]  [--max-chars=max-chars]  [-P max-procs] [--max-procs=max-procs]
       
    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
    转自:https://www.jb51.net/article/103875.htm

    邵先生有话说:先说说我为什么要转这篇文章,因为在之前的工作中每个样本名字都是不一样的,所以不需要特殊处理,而今天遇到了每个不同目录下的样本都是一样的情况,这个处理就需要一些技巧,所以今天用了$(basename $(dirname file)),这样就取出了目录的名字,还有${}这个用法,这些在正常的工作中都有可能会遇到,所以提前预习做好工作很重要!!

    --------------------------------------------------------------------------------------------------------------很多时候在使用Linux的shell时,我们都需要对文件名或目录名进行处理,通常的操作是由路径中提取出文件名,从路径中提取出目录名,提取文件后缀名等等。例如,从路径/dir1/dir2/file.txt中提取也文件名file.txt,提取出目录/dir1/dir2,提取出文件后缀txt等。

    下面介绍两种常用的方法来进行相关的操作。

    一、使用${}

    1、${va

    阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
      

    新浪BLOG意见反馈留言板 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有