2019-10-10-生信人的linux考试后10题-解题记录

11.安装 samtools 软件

生信软件安装：参考https://www.jianshu.com/p/a2a7510908cf

  • 登录服务器
ssh gz16@118.24.223.116
pd19162

1 安装Miniconda3

mkdir -p ~/biosoft
cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
# 查看conda是否按照成功
(base) gz16@VM-0-17-ubuntu:~/biosoft$ conda --version
conda 4.7.10

2.设置conda镜像

source ~/.bashrc
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

3.使用conda 安装samtools 软件

#如果不确定之前是否安装过，可先查看
conda search samtools
#安装指令
conda install -y samtools
#安装过程（省略中间部分）
(base) gz16@VM-0-17-ubuntu:~/biosoft$ conda install -y samtools
Collecting package metadata (current_repodata.json): done
Solving environment: done
...
Downloading and Extracting Packages
libssh2-1.8.2        | 257 KB    | ##################################### | 100%
samtools-1.9         | 299 KB    | ##################################### | 100%
certifi-2019.9.11    | 147 KB    | ##################################### | 100%
ca-certificates-2019 | 144 KB    | ##################################### | 100%
openssl-1.1.1c       | 2.1 MB    | ##################################### | 100%
conda-4.7.12         | 3.0 MB    | ##################################### | 100%
htslib-1.9           | 1.2 MB    | ##################################### | 100%
libdeflate-1.2       | 63 KB     | ##################################### | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

十二、打开 后缀为BAM 的文件，找到产生该文件的命令。 提示一下命令是：

/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2-align-s --wrapper basic-0 -p 20 -x /home/jianmingzeng/reference/index/bowtie/hg38 -S /home/jianmingzeng/data/public/allMouse/alignment/WT_rep2_Input.sam -U /tmp/41440.unp

由于在当前服务器没有找到bam文件，查到前10题里面有个下载的压缩包里面有bam文件
返回复现

九、下载http://www.biotrainee.com/jmzeng/rmDuplicate.zip文件，并且解压，查看里面的文件夹结构。

wget http://www.biotrainee.com/jmzeng/rmDuplicate.zip
ls -ll
unzip rmDuplicate.zip
cd rmDuplicate
tree

十、打开第九题解压的文件，进入rmDuplicate/samtools/single文件夹里面，查看后缀为.sam的文件，搞清楚 生物信息学里面的SAM/BAM定义是什么。

cd samtools/single/
less -S tmp.sam #最佳查看方式
cat tmp.sam
vim tmp.sam
#vim编辑器退出需要先`esc`，输入`:`输入`q/q!/wq`退出
# 打开bam 文件前2行
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ less -S tmp.sam | head -n 2
SRR1042600.42157053 0   chr1    629895  42  51M *   0   0   ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA CCCFFFFFHHHHHJJJJJJJ#4AGHJJIIJJIIIIIJJJJIJIIIIJJIJI AS:i:-6 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11C8A30    YT:Z:UU
SRR1042600.42212881 0   chr1    629895  42  51M *   0   0   ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA @@<FDFFBFDHHFJEIIGJI#3AFHGEHEIJIIGIIGGIJIIJIGIIGIIJ AS:i:-6 XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:11C8A30    YT:Z:UU

十三题、根据上面的命令，找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38具体有多少条染色体。

由于当前服务器没有/home/jianmingzeng/reference/index/bowtie/hg38所有还是在上面的题目中解压的文件目录下查找
找到在~/rmDuplicate/samtools/single目录下游2个bam文件

(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ ls
readme.txt  tmp.rmdup.bam     tmp.sam         tmp.sorted.vcf.gz
tmp.header  tmp.rmdup.vcf.gz  tmp.sorted.bam
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ samtools view -H tmp.rmdup.bam | head -n 10
@HD VN:1.0  SO:coordinate
@SQ SN:chr1 LN:248956422
@SQ SN:chr10    LN:133797422
@SQ SN:chr11    LN:135086622
@SQ SN:chr11_KI270721v1_random  LN:100316
@SQ SN:chr12    LN:133275309
@SQ SN:chr13    LN:114364328
@SQ SN:chr14    LN:107043718
@SQ SN:chr14_GL000009v2_random  LN:201709
@SQ SN:chr14_GL000225v1_random  LN:211173
#因为是二进制压缩文件，查看不能用less而应该用zless
zless -S tmp.rmdup.bam

查看生信基础数据格式参考生信菜鸟团-专题学习目录(2)
数据格式专题
1.FastQ和FastA格式
2.SAM格式
3.gff/gtf格式
4.Bigwig/Wiggle格式
5.bed格式
6.vcf格式
7.Blast&Blat

简单探索后找到染色体信息

(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ samtools view -H tmp.rmdup.bam | awk '{print $2}'| sort | uniq -c| grep -v '_'
      1 ID:bowtie2
      1 SN:chr1
       ...
      1 SN:chrM
      1 SN:chrX
      1 SN:chrY
      1 VN:1.0

十四题、上面的后缀为bam的文件的第二列，只有 0 和 16 两个数字，用 cut/sort/uniq等命令统计它们的个数。

因为和原来的数据不同，这里tmp.rmdup.bam文件第二列是LN:248956422应该是染色体位置信息但不是单纯的数值，因此换别的bam文件

(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ cat tmp.sam | awk '{print $2}' | sort | uniq -c
     29 0
     24 16
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/single$ samtools view tmp.sorted.bam |awk '{print $2}' | sort | uniq -c
     29 0
     24 16

十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件，再次查看第二列，并且统计。

(base) gz16@VM-0-17-ubuntu:~$ cd rmDuplicate/samtools/paired/
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ ls
readme.txt  tmp.rmdup.bam     tmp.sam         tmp.sorted.vcf.gz
tmp.header  tmp.rmdup.vcf.gz  tmp.sorted.bam
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ zless -S tmp.sorted.bam | head -n 2
BAM'
    @HD VN:1.3  SO:coordinate
@SQ SN:chr10    LN:135534747
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ samtools view tmp.sorted.bam | head -n 2
D00691:39:C7HGRANXX:7:1102:7445:18770   99  chr10   93614   60  126M    =93621  133 GGCACGTGGTGACCCCACTCATGGTAGCAGACACCAGGTGGTTCAGGTCACCATAGGTGGGTGTGGGCAGTTTTAGGGTCTTGGAACATATGTCATACAGAGCTTCGTTATCTATGCAAAAGGTCT  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF<FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NM:i:0  MD:Z:126    AS:i:126    XS:i:111XA:Z:chr3,-197846908,126M,3;chr9,-141070974,126M,3;chr1,-808987,126M,4;chr4,+190904265,126M,4;  MQ:i:60
D00691:39:C7HGRANXX:7:1102:7445:18770   147 chr10   93621   60  126M    =93614  -133    GGTGACCCCACTCATGGTAGCAGACACCAGGTGGTTCAGGTCACCATAGGTGGGTGTGGGCAGTTTTAGGGTCTTGGAACATATGTCATACAGAGCTTCGTTATCTATGCAAAAGGTCTCATCTGC  FFFFFFFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFFFFFFFFFFBBBBB  NM:i:0  MD:Z:126    AS:i:126    XS:i:111XA:Z:chr9,+141070967,126M,3;chr3,+197846902,1S125M,3;   MQ:i:60
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ samtools view tmp.sorted.bam | awk '{print $2}'
99
147
323
387
353
97
371
433
97
99
147
99
147
99
147
99
147
99
99
147
147
163
163
83
83
163
83
99
147
99
(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ samtools view tmp.sorted.bam | awk '{print $2}'| sort |uniq -c
      8 147
      3 163
      1 323
      1 353
      1 371
      1 387
      1 433
      3 83
      2 97
      9 99

十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，查看里面的文件夹结构， 这个文件有2.3M，注意留心下载时间及下载速度。

(base) gz16@VM-0-17-ubuntu:~/rmDuplicate/samtools/paired$ cd ~
(base) gz16@VM-0-17-ubuntu:~$ wget http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
--2019-10-10 12:19:16--  http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
Resolving www.biotrainee.com (www.biotrainee.com)... 123.206.72.184
Connecting to www.biotrainee.com (www.biotrainee.com)|123.206.72.184|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 2391084 (2.3M) [application/zip]
Saving to: 'sickle-results.zip'

sickle-results.zip   32%[=====>              ] 765.01K  51.3KB/s    eta 21s

#下载完成后解压缩并查看目录结构
2019-10-10 12:19:36 (118 KB/s) - 'sickle-results.zip' saved [2391084/2391084]

(base) gz16@VM-0-17-ubuntu:~$ unzip sickle-results.zip
Archive:  sickle-results.zip
   creating: sickle-results/
  inflating: sickle-results/command.txt
  ...
  inflating: sickle-results/trimmed_output_file2_fastqc.zip
(base) gz16@VM-0-17-ubuntu:~$ cd sickle-results
(base) gz16@VM-0-17-ubuntu:~/sickle-results$ tree
.
|-- command.txt
|-- single_tmp_fastqc.html
|-- single_tmp_fastqc.zip
|-- test1_fastqc.html
|-- test1_fastqc.zip
|-- test2_fastqc.html
|-- test2_fastqc.zip
|-- trimmed_output_file1_fastqc.html
|-- trimmed_output_file1_fastqc.zip
|-- trimmed_output_file2_fastqc.html
`-- trimmed_output_file2_fastqc.zip

这里面有个command.txt,其中有些代码可参考安装软件

十七题、解压sickle-results/single_tmp_fastqc.zip文件，并且进入解压后的文件夹，找到fastqc_data.txt文件，并且搜索该文本文件以 >>开头的有多少行？

插曲，忘记基础解压缩指令，反而写成压缩指令，回头复习一下

  • 1).gzip 和gunzip 指令
    功能： gzip 用于压缩文件
    gunzip 用于解压文件
  • 2). zip指令 和 unzip 指令
    功能：zip 用于压缩文件
    unzip 用于解压文件
  • zip 常用选项：
    -r 递归压缩
    unzip 常用选项：
    -d <目录> 指定解压后的文件的存放目录
  • 3). tar 指令
    功能：打包指令最后打包的文件是.tar.gz文件
    注意：打包用-zcvf
    解压用-zxvf
(base) gz16@VM-0-17-ubuntu:~/sickle-results$ gunzip single_tmp_fastqc.zip.gz
(base) gz16@VM-0-17-ubuntu:~/sickle-results$ ll
total 3028
drwxrwxr-x  2 gz16 gz16   4096 Oct 10 12:33 ./
drwxr-xr-x 14 gz16 gz     4096 Oct 10 12:20 ../
-rw-rw-r--  1 gz16 gz16    901 Oct  6  2016 command.txt
-rw-rw-r--  1 gz16 gz16 259290 Oct  6  2016 single_tmp_fastqc.html
-rw-rw-r--  1 gz16 gz16 260908 Oct  6  2016 single_tmp_fastqc.zip
-rw-rw-r--  1 gz16 gz16 308721 Oct  6  2016 test1_fastqc.html
-rw-rw-r--  1 gz16 gz16 334674 Oct  6  2016 test1_fastqc.zip
-rw-rw-r--  1 gz16 gz16 305348 Oct  6  2016 test2_fastqc.html
-rw-rw-r--  1 gz16 gz16 332699 Oct  6  2016 test2_fastqc.zip
-rw-rw-r--  1 gz16 gz16 305987 Oct  6  2016 trimmed_output_file1_fastqc.html
-rw-rw-r--  1 gz16 gz16 329011 Oct  6  2016 trimmed_output_file1_fastqc.zip
-rw-rw-r--  1 gz16 gz16 302503 Oct  6  2016 trimmed_output_file2_fastqc.html
-rw-rw-r--  1 gz16 gz16 328552 Oct  6  2016 trimmed_output_file2_fastqc.zip
(base) gz16@VM-0-17-ubuntu:~/sickle-results$ unzip single_tmp_fastqc.zip
Archive:  single_tmp_fastqc.zip
   creating: single_tmp_fastqc/
   creating: single_tmp_fastqc/Icons/
   creating: single_tmp_fastqc/Images/
  inflating: single_tmp_fastqc/Icons/fastqc_icon.png
  ...
  inflating: single_tmp_fastqc/fastqc.fo
(base) gz16@VM-0-17-ubuntu:~/sickle-results$ cd single_tmp_fastqc/
(base) gz16@VM-0-17-ubuntu:~/sickle-results/single_tmp_fastqc$ ls
Icons  Images  fastqc.fo  fastqc_data.txt  fastqc_report.html  summary.txt
(base) gz16@VM-0-17-ubuntu:~/sickle-results/single_tmp_fastqc$ less -S fastqc_data.txt | awk '/^>>/{print $0}'
>>Basic Statistics  pass
>>END_MODULE
...
>>Kmer Content  warn
>>END_MODULE
(base) gz16@VM-0-17-ubuntu:~/sickle-results/single_tmp_fastqc$ less -S fastqc_data.txt | awk '/^>>/{print $0}' | wc -l
24

十八题、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件，去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID，然后找到它们的hg38.tss 文件的哪一行。

https://www.ncbi.nlm.nih.gov/gene/1827

cd ~
wget http://www.biotrainee.com/jmzeng/tmp/hg38.tss

(base) gz16@VM-0-17-ubuntu:~$ head hg38.tss
NR_046018   chr1    9874    13874   0
NR_024540   chr1    27370   31370   1
NR_104148   chr7    64664083    64668083    0
NR_111960   chrX    44871175    44875175    0
NR_028458   chr14   92104621    92108621    1
NR_028459   chr14   92104621    92108621    1
NR_026818   chr1    34081   38081   1
NR_026820   chr1    34081   38081   1
NR_026822   chr1    34081   38081   1
NM_001005484    chr1    67091   71091   0
#参考基因NR_...或者NM_...
#在NCBI选取目的基因RCAN1转录本c:NM_203418 在hg38.tss文件的第3467行
(base) gz16@VM-0-17-ubuntu:~$ grep -n NM_203418 hg38.tss
3467:NM_203418  chr21   34525010    34529010    1

太棒了，原来还有这种操作，我发现又有一个问题解决了，我需要找到6号染色体体短臂的一段序列有哪些基因，这就可以提取出来基因名称了呀，棒棒的！

十九题、解析hg38.tss文件，统计每条染色体的基因个数。

(base) gz16@VM-0-17-ubuntu:~$ cat hg38.tss|awk '{print $2}' |sort |uniq -c
   6050 chr1
   2824 chr10
      2 chr10_GL383545v1_alt
     10 chr10_GL383546v1_alt
      2 chr10_KI270825v1_alt
   3449 chr11
  ...
   2553 chrX
      3 chrX_KI270880v1_alt
      4 chrX_KI270881v1_alt
      1 chrX_KI270913v1_alt
    414 chrY

这样来看应该是染色体还有很多的变体（有些不是完整的染色体），这里就需要在去查查资料关于染色体的相关背景知识了，解答此题需要先取出这些不完整的染色体

(base) gz16@VM-0-17-ubuntu:~$ cat hg38.tss|awk '{print $2}' | sort | uniq -c | grep -v '_'
   6050 chr1
   2824 chr10
   3449 chr11
   2931 chr12
   1122 chr13
   1883 chr14
   2168 chr15
   2507 chr16
   3309 chr17
    873 chr18
   3817 chr19
   4042 chr2
   1676 chr20
    868 chr21
   1274 chr22
   3277 chr3
   2250 chr4
   2684 chr5
   3029 chr6
   2720 chr7
   2069 chr8
   2301 chr9
      2 chrM
   2553 chrX
    414 chrY

发现还有2 个基因在chrM上，这个是什么染色体，需要补一下背景知识

二十题、解析hg38.tss文件，统计NM和NR开头的数量，了解NM和NR开头的含义

(base) gz16@VM-0-17-ubuntu:~$ cat hg38.tss |awk '{print $1}'|cut -c 1-2 |sort| uniq -c
  51064 NM
  15954 NR

关于NM和NR开头的含义还是需要查看生信菜鸟团的基础知识专题。

推荐
生信技能树公益视频合辑
生信技能树简书账号
生信工程师入门最佳指南
生信技能树全球公益巡讲
招学徒

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 159,835评论 4赞 364
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 67,598评论 1赞 295
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 109,569评论 0赞 244
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 44,159评论 0赞 213
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 52,533评论 3赞 287
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 40,710评论 1赞 222
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 31,923评论 2赞 313
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 30,674评论 0赞 203
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 34,421评论 1赞 246
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,622评论 2赞 245
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,115评论 1赞 260
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 28,428评论 2赞 254
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 33,114评论 3赞 238
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,097评论 0赞 8
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,875评论 0赞 197
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 35,753评论 2赞 276
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 35,649评论 2赞 271