生信入门之初级作业及答案

2023年3月25日 229次阅读 来源: 秦城听雪

1. 在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列。

mkdir -p 1/2/3/4/5/6/7/8/9

2. 在创建好的文件夹下面，比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ，里面创建文本文件 me.txt

touch /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9/me.txt

3. 在文本文件 me.txt 里面输入内容:

Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?

方法一：

vi me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
:wq

方法二：

cat> me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
##按ctrl + d保存退出

4.删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt

rm -r 1/2/3/4/5/6/7/8/9/me.txt

5.在任意文件夹下面创建 folder1~5这5个文件夹，然后每个文件夹下面继续创建 folder1~5这5个文件夹

mkdir -p  folder{1..5}/folder{1..5}

6.在第五题创建的每一个文件夹下面都创建第二题文本文件 me.txt ，内容也要一样。

cat >me.txt 
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
#按ctrl + d保存退出

for a in {1..5};do
     cd ~/folder$a
            for b in {1..5};do
                    cd ~/folder$a/folder$b
  cp ./me.txt  ~/folder$a/folder$b
     done
done

7.再次删除掉前面几个步骤建立的文件夹及文件

rm -r folder{1..5}

8. 下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件，后在里面选择含有 H3K4me3 的那一行是第几行，该文件总共有几行。

wget  -c http://www.biotrainee.com/jmzeng/igv/test.bed
grep -n  'H3K4me3'  test.bed
wc -l test.bed

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

wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip
unzip rmDuplicate.zip
tree rmDuplicate

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

cd rmDuplicate/samtools/single
ls *.sam

SAM(Sequence Alignment/Map)格式是一种通用的比对格式，用来存储reads到参考序列的比对信息。
BAM是SAM的二进制格式，因此两者格式相同，只是BAM文件占用储存空间更小，运算更快。

11.安装 samtools 软件

利用conda安装

source ~/miniconda3/bin/activate
conda search samtools
conda install samtools

12.打开后缀为BAM 的文件，找到产生该文件的命令。

#查看一个文件名叫做SRR1039510.sort.bam文件的产生命令
samtools view -h SRR1039510.sort.bam |grep '^@PG'|awk 'BEGIN{FS="\t"}{print $5}'|cut -d: -f2

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

#查看   注：此处用的文件为SRR1039510.sort.bam文件
samtools view -h SRR1039510.sort.bam |egrep '^@S.*?(chr[XYM]\s+.*|chr[1-9]?[0-9]\s+).*'|less
#计数
samtools view -h SRR1039510.sort.bam |egrep '^@S.*?(chr[XYM]\s+.*|chr[1-9]?[0-9]\s+).*'|wc -l

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

samtools view SRR1039510.sort.bam |cut -f2|sort |uniq -c

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

samtools view SRR1039510.sort.bam|cut -f2 |sort |uniq -c|sort -t' ' -nrk1,1

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

wget http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
unzip sickle-results
tree

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

cd sickle-results
unzip single_tmp_fastqc.zip
cd single_tmp_fastqc
search  '^>>'  fastqc_data.txt |wc -l

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

wget http://www.biotrainee.com/jmzeng/tmp/hg38.tss
grep 'NM_000546' hg38.tss 
grep 'NM_001126113' hg38.tss

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

cat hg38.tss |cut -f2|sort|uniq -c

20.解析hg38.tss 文件，统计NM和NR开头的序列，了解NM和NR开头的含义。

grep '^NM' hg38.tss |wc -l
grep '^NR' hg38.tss |wc -l

或

grep -oE '^(NM|NR)' hg38.tss |sort|uniq -c

NM：转录组产物的序列mRNA
NR：非编码的转录组序列ncRNA

    原文作者：秦城听雪
    原文地址: https://www.jianshu.com/p/da43ef737f74
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。
点赞